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

Improve I/O on TC integrals

Browse Source
This commit is contained in:
Anthony Scemama 2023-03-17 19:23:07 +01:00
parent 29230b175d
commit e710d26237
5 changed files with 118 additions and 193 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
src/non_h_ints_mu/grad_squared.irp.f
View File

@ -2,7 +2,7 @@
! --- ! ---
! TODO : strong optmization : write the loops in a different way ! TODO : strong optmization : write the loops in a different way
! : for each couple of AO, the gaussian product are done once for all ! : for each couple of AO, the gaussian product are done once for all
BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_points_final_grid) ] BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_points_final_grid) ]
@ -20,14 +20,14 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi
! gradu_squared_u_ij_mu = -0.50 x \int r2 \phi_i(2) \phi_j(2) [ v1^2 v2^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) + u12^2 v2^2 (grad_1 v1)^2 + 2 u12 v1 v2^2 (grad_1 u12) . (grad_1 v1) ] ! gradu_squared_u_ij_mu = -0.50 x \int r2 \phi_i(2) \phi_j(2) [ v1^2 v2^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) + u12^2 v2^2 (grad_1 v1)^2 + 2 u12 v1 v2^2 (grad_1 u12) . (grad_1 v1) ]
! = -0.25 x v1^2 \int r2 \phi_i(2) \phi_j(2) [1 - erf(mu r12)]^2 v2^2 ! = -0.25 x v1^2 \int r2 \phi_i(2) \phi_j(2) [1 - erf(mu r12)]^2 v2^2
! + -0.50 x (grad_1 v1)^2 \int r2 \phi_i(2) \phi_j(2) u12^2 v2^2 ! + -0.50 x (grad_1 v1)^2 \int r2 \phi_i(2) \phi_j(2) u12^2 v2^2
! + -1.00 x v1 (grad_1 v1) \int r2 \phi_i(2) \phi_j(2) (grad_1 u12) v2^2 ! + -1.00 x v1 (grad_1 v1) \int r2 \phi_i(2) \phi_j(2) (grad_1 u12) v2^2
! = v1^2 x int2_grad1u2_grad2u2_j1b2 ! = v1^2 x int2_grad1u2_grad2u2_j1b2
! + -0.5 x (grad_1 v1)^2 x int2_u2_j1b2 ! + -0.5 x (grad_1 v1)^2 x int2_u2_j1b2
! + -1.0 X V1 x (grad_1 v1) \cdot [ int2_u_grad1u_j1b2 x r - int2_u_grad1u_x_j1b ] ! + -1.0 X V1 x (grad_1 v1) \cdot [ int2_u_grad1u_j1b2 x r - int2_u_grad1u_x_j1b ]
! !
! !
END_DOC END_DOC
implicit none implicit none
integer :: ipoint, i, j, m, igauss integer :: ipoint, i, j, m, igauss
double precision :: x, y, z, r(3), delta, coef double precision :: x, y, z, r(3), delta, coef
@ -100,7 +100,7 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for gradu_squared_u_ij_mu = ', time1 - time0 print*, ' Wall time for gradu_squared_u_ij_mu = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---
@ -151,7 +151,7 @@ END_PROVIDER
! !
! deallocate(ac_mat) ! deallocate(ac_mat)
! !
!END_PROVIDER !END_PROVIDER
! --- ! ---
@ -214,12 +214,12 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_loop, (ao_num, ao_num, ao_nu
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for tc_grad_square_ao_loop = ', time1 - time0 print*, ' Wall time for tc_grad_square_ao_loop = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_grid) ] BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_grid) ]
implicit none implicit none
integer :: ipoint, i, j, m, igauss integer :: ipoint, i, j, m, igauss
double precision :: r(3), delta, coef double precision :: r(3), delta, coef
@ -267,7 +267,7 @@ BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_g
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for grad12_j12 = ', time1 - time0 print*, ' Wall time for grad12_j12 = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---
@ -297,12 +297,12 @@ BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for u12sq_j1bsq = ', time1 - time0 print*, ' Wall time for u12sq_j1bsq = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num, n_points_final_grid) ] BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num, n_points_final_grid) ]
implicit none implicit none
integer :: ipoint, i, j, m, igauss integer :: ipoint, i, j, m, igauss
double precision :: x, y, z double precision :: x, y, z
@ -347,7 +347,7 @@ BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num,
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b = ', time1 - time0 print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---
@ -370,26 +370,18 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
if(read_tc_integ) then if(read_tc_integ) then
open(unit=11, form="unformatted", file='tc_grad_square_ao', action="read") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_square_ao', action="read")
do i = 1, ao_num read(11) tc_grad_square_ao
do j = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
read(11) tc_grad_square_ao(l,k,j,i)
enddo
enddo
enddo
enddo
close(11) close(11)
else else
allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid)) allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid))
b_mat = 0.d0 b_mat = 0.d0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, k, ipoint) & !$OMP PRIVATE (i, k, ipoint) &
!$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector) !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do i = 1, ao_num do i = 1, ao_num
@ -401,11 +393,11 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
tmp = 0.d0 tmp = 0.d0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j, l, ipoint) & !$OMP PRIVATE (j, l, ipoint) &
!$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12) !$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -417,25 +409,25 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
tc_grad_square_ao = 0.d0 tc_grad_square_ao = 0.d0
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
, tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & , tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
, 1.d0, tc_grad_square_ao, ao_num*ao_num) , 1.d0, tc_grad_square_ao, ao_num*ao_num)
deallocate(tmp, b_mat) deallocate(tmp, b_mat)
call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num) call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num)
!!$OMP PARALLEL & !!$OMP PARALLEL &
!!$OMP DEFAULT (NONE) & !!$OMP DEFAULT (NONE) &
!!$OMP PRIVATE (i, j, k, l) & !!$OMP PRIVATE (i, j, k, l) &
!!$OMP SHARED (ac_mat, tc_grad_square_ao, ao_num) !!$OMP SHARED (ac_mat, tc_grad_square_ao, ao_num)
!!$OMP DO SCHEDULE (static) !!$OMP DO SCHEDULE (static)
! do j = 1, ao_num ! do j = 1, ao_num
! do l = 1, ao_num ! do l = 1, ao_num
! do i = 1, ao_num ! do i = 1, ao_num
! do k = 1, ao_num ! do k = 1, ao_num
! tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) ! tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
! enddo ! enddo
! enddo ! enddo
! enddo ! enddo
@ -444,23 +436,17 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
!!$OMP END PARALLEL !!$OMP END PARALLEL
endif endif
if(write_tc_integ) then if(write_tc_integ.and.mpi_master) then
open(unit=11, form="unformatted", file='tc_grad_square_ao', action="write") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_square_ao', action="write")
do i = 1, ao_num call ezfio_set_work_empty(.False.)
do j = 1, ao_num write(11) tc_grad_square_ao
do k = 1, ao_num
do l = 1, ao_num
write(11) tc_grad_square_ao(l,k,j,i)
enddo
enddo
enddo
enddo
close(11) close(11)
call ezfio_set_tc_keywords_io_tc_integ('Read')
endif endif
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for tc_grad_square_ao = ', time1 - time0 print*, ' Wall time for tc_grad_square_ao = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---

76
src/non_h_ints_mu/grad_squared_manu.irp.f
View File

@ -17,29 +17,21 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
call wall_time(time0) call wall_time(time0)
if(read_tc_integ) then if(read_tc_integ) then
open(unit=11, form="unformatted", file='tc_grad_square_ao_test', action="read") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_square_ao_test', action="read")
do i = 1, ao_num read(11) tc_grad_square_ao_test
do j = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
read(11) tc_grad_square_ao_test(l,k,j,i)
enddo
enddo
enddo
enddo
close(11) close(11)
else else
provide u12sq_j1bsq_test u12_grad1_u12_j1b_grad1_j1b_test grad12_j12_test provide u12sq_j1bsq_test u12_grad1_u12_j1b_grad1_j1b_test grad12_j12_test
allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid)) allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid))
b_mat = 0.d0 b_mat = 0.d0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, k, ipoint) & !$OMP PRIVATE (i, k, ipoint) &
!$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector) !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do i = 1, ao_num do i = 1, ao_num
@ -51,11 +43,11 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
tmp = 0.d0 tmp = 0.d0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j, l, ipoint) & !$OMP PRIVATE (j, l, ipoint) &
!$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq_test, u12_grad1_u12_j1b_grad1_j1b_test, grad12_j12_test) !$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq_test, u12_grad1_u12_j1b_grad1_j1b_test, grad12_j12_test)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -67,23 +59,23 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
tc_grad_square_ao_test = 0.d0 tc_grad_square_ao_test = 0.d0
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
, tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & , tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
, 1.d0, tc_grad_square_ao_test, ao_num*ao_num) , 1.d0, tc_grad_square_ao_test, ao_num*ao_num)
deallocate(tmp, b_mat) deallocate(tmp, b_mat)
call sum_A_At(tc_grad_square_ao_test(1,1,1,1), ao_num*ao_num) call sum_A_At(tc_grad_square_ao_test(1,1,1,1), ao_num*ao_num)
!do i = 1, ao_num !do i = 1, ao_num
! do j = 1, ao_num ! do j = 1, ao_num
! do k = i, ao_num ! do k = i, ao_num
! do l = max(j,k), ao_num ! do l = max(j,k), ao_num
! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j)) ! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j))
! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l) ! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l)
! end do ! end do
! !if (j.eq.k) then ! !if (j.eq.k) then
! ! do l = j+1, ao_num ! ! do l = j+1, ao_num
! ! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j)) ! ! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j))
@ -95,14 +87,14 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
! ! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l) ! ! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l)
! ! enddo ! ! enddo
! !endif ! !endif
! enddo ! enddo
! enddo ! enddo
!enddo !enddo
!tc_grad_square_ao_test = 2.d0 * tc_grad_square_ao_test !tc_grad_square_ao_test = 2.d0 * tc_grad_square_ao_test
! !$OMP PARALLEL & ! !$OMP PARALLEL &
! !$OMP DEFAULT (NONE) & ! !$OMP DEFAULT (NONE) &
! !$OMP PRIVATE (i, j, k, l) & ! !$OMP PRIVATE (i, j, k, l) &
! !$OMP SHARED (tc_grad_square_ao_test, ao_num) ! !$OMP SHARED (tc_grad_square_ao_test, ao_num)
! !$OMP DO SCHEDULE (static) ! !$OMP DO SCHEDULE (static)
! integer :: ii ! integer :: ii
@ -121,10 +113,10 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
! print *, ' ii =', ii ! print *, ' ii =', ii
! !$OMP END DO ! !$OMP END DO
! !$OMP END PARALLEL ! !$OMP END PARALLEL
! !$OMP PARALLEL & ! !$OMP PARALLEL &
! !$OMP DEFAULT (NONE) & ! !$OMP DEFAULT (NONE) &
! !$OMP PRIVATE (i, j, k, l) & ! !$OMP PRIVATE (i, j, k, l) &
! !$OMP SHARED (tc_grad_square_ao_test, ao_num) ! !$OMP SHARED (tc_grad_square_ao_test, ao_num)
! !$OMP DO SCHEDULE (static) ! !$OMP DO SCHEDULE (static)
! do j = 1, ao_num ! do j = 1, ao_num
@ -144,24 +136,18 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
endif endif
if(write_tc_integ) then if(write_tc_integ.and.mpi_master) then
open(unit=11, form="unformatted", file='tc_grad_square_ao_test', action="write") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_square_ao_test', action="write")
do i = 1, ao_num call ezfio_set_work_empty(.False.)
do j = 1, ao_num write(11) tc_grad_square_ao_test
do k = 1, ao_num
do l = 1, ao_num
write(11) tc_grad_square_ao_test(l,k,j,i)
enddo
enddo
enddo
enddo
close(11) close(11)
call ezfio_set_tc_keywords_io_tc_integ('Read')
endif endif
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for tc_grad_square_ao_test = ', time1 - time0 print*, ' Wall time for tc_grad_square_ao_test = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---
@ -189,7 +175,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test_ref, (ao_num, ao_num, a
b_mat = 0.d0 b_mat = 0.d0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, k, ipoint) & !$OMP PRIVATE (i, k, ipoint) &
!$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector) !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do i = 1, ao_num do i = 1, ao_num
@ -205,7 +191,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test_ref, (ao_num, ao_num, a
tmp = 0.d0 tmp = 0.d0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j, l, ipoint) & !$OMP PRIVATE (j, l, ipoint) &
!$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq_test, u12_grad1_u12_j1b_grad1_j1b_test, grad12_j12_test) !$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq_test, u12_grad1_u12_j1b_grad1_j1b_test, grad12_j12_test)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -226,7 +212,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test_ref, (ao_num, ao_num, a
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l) & !$OMP PRIVATE (i, j, k, l) &
!$OMP SHARED (ac_mat, tc_grad_square_ao_test_ref, ao_num) !$OMP SHARED (ac_mat, tc_grad_square_ao_test_ref, ao_num)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do j = 1, ao_num do j = 1, ao_num
@ -246,7 +232,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test_ref, (ao_num, ao_num, a
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for tc_grad_square_ao_test_ref = ', time1 - time0 print*, ' Wall time for tc_grad_square_ao_test_ref = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---
@ -276,12 +262,12 @@ BEGIN_PROVIDER [ double precision, u12sq_j1bsq_test, (ao_num, ao_num, n_points_f
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for u12sq_j1bsq_test = ', time1 - time0 print*, ' Wall time for u12sq_j1bsq_test = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b_test, (ao_num, ao_num, n_points_final_grid) ] BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b_test, (ao_num, ao_num, n_points_final_grid) ]
implicit none implicit none
integer :: ipoint, i, j, m, igauss integer :: ipoint, i, j, m, igauss
double precision :: x, y, z double precision :: x, y, z
@ -328,12 +314,12 @@ BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b_test, (ao_num, ao
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b_test = ', time1 - time0 print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b_test = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, grad12_j12_test, (ao_num, ao_num, n_points_final_grid) ] BEGIN_PROVIDER [ double precision, grad12_j12_test, (ao_num, ao_num, n_points_final_grid) ]
implicit none implicit none
integer :: ipoint, i, j, m, igauss integer :: ipoint, i, j, m, igauss
double precision :: r(3), delta, coef double precision :: r(3), delta, coef
@ -381,7 +367,7 @@ BEGIN_PROVIDER [ double precision, grad12_j12_test, (ao_num, ao_num, n_points_fi
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for grad12_j12_test = ', time1 - time0 print*, ' Wall time for grad12_j12_test = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---

28
src/non_h_ints_mu/new_grad_tc.irp.f
View File

@ -36,16 +36,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_
if(read_tc_integ) then if(read_tc_integ) then
open(unit=11, form="unformatted", file='int2_grad1_u12_ao', action="read") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read")
do m = 1, 3 read(11) int2_grad1_u12_ao
do ipoint = 1, n_points_final_grid
do j = 1, ao_num
do i = 1, ao_num
read(11) int2_grad1_u12_ao(i,j,ipoint,m)
enddo
enddo
enddo
enddo
close(11) close(11)
else else
@ -89,18 +81,12 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_
endif endif
if(write_tc_integ) then if(write_tc_integ.and.mpi_master) then
open(unit=11, form="unformatted", file='int2_grad1_u12_ao', action="write") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="write")
do m = 1, 3 call ezfio_set_work_empty(.False.)
do ipoint = 1, n_points_final_grid write(11) int2_grad1_u12_ao
do j = 1, ao_num
do i = 1, ao_num
write(11) int2_grad1_u12_ao(i,j,ipoint,m)
enddo
enddo
enddo
enddo
close(11) close(11)
call ezfio_set_tc_keywords_io_tc_integ('Read')
endif endif
call wall_time(time1) call wall_time(time1)

103
src/non_h_ints_mu/new_grad_tc_manu.irp.f
View File

@ -3,7 +3,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
BEGIN_DOC BEGIN_DOC
! !
! int2_grad1_u12_ao_test(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) ! int2_grad1_u12_ao_test(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2)
! !
! where r1 = r(ipoint) ! where r1 = r(ipoint)
! !
@ -15,9 +15,9 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
! if J(r1,r2) = u12 x v1 x v2 ! if J(r1,r2) = u12 x v1 x v2
! !
! int2_grad1_u12_ao_test(i,j,ipoint,:) = v1 x [ 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] v2 \phi_i(r2) \phi_j(r2) ] ! int2_grad1_u12_ao_test(i,j,ipoint,:) = v1 x [ 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] v2 \phi_i(r2) \phi_j(r2) ]
! - \grad_1 v1 x [ \int dr2 u12 v2 \phi_i(r2) \phi_j(r2) ] ! - \grad_1 v1 x [ \int dr2 u12 v2 \phi_i(r2) \phi_j(r2) ]
! = 0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) ! = 0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:)
! - 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) ! - 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:)
! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint) ! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
! !
! !
@ -35,25 +35,18 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
if(read_tc_integ) then if(read_tc_integ) then
open(unit=11, form="unformatted", file='int2_grad1_u12_ao_test', action="read") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao_test', action="read")
do m = 1, 3 read(11) int2_grad1_u12_ao_test
do ipoint = 1, n_points_final_grid
do j = 1, ao_num
do i = 1, ao_num
read(11) int2_grad1_u12_ao_test(i,j,ipoint,m)
enddo
enddo
enddo
enddo
close(11) close(11)
else else
if(j1b_type .eq. 3) then if(j1b_type .eq. 3) then
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint) x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint) y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint) z = final_grid_points(3,ipoint)
tmp0 = 0.5d0 * v_1b(ipoint) tmp0 = 0.5d0 * v_1b(ipoint)
tmp_x = v_1b_grad(1,ipoint) tmp_x = v_1b_grad(1,ipoint)
tmp_y = v_1b_grad(2,ipoint) tmp_y = v_1b_grad(2,ipoint)
@ -87,24 +80,18 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
endif endif
if(write_tc_integ) then if(write_tc_integ.and.mpi_master) then
open(unit=11, form="unformatted", file='int2_grad1_u12_ao_test', action="write") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao_test', action="write")
do m = 1, 3 call ezfio_set_work_empty(.False.)
do ipoint = 1, n_points_final_grid write(11) int2_grad1_u12_ao_test
do j = 1, ao_num
do i = 1, ao_num
write(11) int2_grad1_u12_ao_test(i,j,ipoint,m)
enddo
enddo
enddo
enddo
close(11) close(11)
call ezfio_set_tc_keywords_io_tc_integ('Read')
endif endif
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for int2_grad1_u12_ao_test = ', time1 - time0 print*, ' Wall time for int2_grad1_u12_ao_test = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---
@ -114,9 +101,9 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
! !
! tc_grad_and_lapl_ao_test(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) | ij > ! tc_grad_and_lapl_ao_test(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) | ij >
! !
! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2)
! !
! This is obtained by integration by parts. ! This is obtained by integration by parts.
! !
END_DOC END_DOC
@ -131,40 +118,32 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
call wall_time(time0) call wall_time(time0)
if(read_tc_integ) then if(read_tc_integ) then
open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao_test', action="read") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_and_lapl_ao_test', action="read")
do i = 1, ao_num read(11) tc_grad_and_lapl_ao_test
do j = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
read(11) tc_grad_and_lapl_ao_test(l,k,j,i)
enddo
enddo
enddo
enddo
close(11) close(11)
else else
provide int2_grad1_u12_ao_test provide int2_grad1_u12_ao_test
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num)) allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num))
b_mat = 0.d0 b_mat = 0.d0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) & !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) &
!$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, & !$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, &
!$OMP ao_num, n_points_final_grid, final_weight_at_r_vector) !$OMP ao_num, n_points_final_grid, final_weight_at_r_vector)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do i = 1, ao_num do i = 1, ao_num
do k = 1, ao_num do k = 1, ao_num
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
ao_i_r = aos_in_r_array_transp(ipoint,i) ao_i_r = aos_in_r_array_transp(ipoint,i)
ao_k_r = aos_in_r_array_transp(ipoint,k) ao_k_r = aos_in_r_array_transp(ipoint,k)
b_mat(ipoint,k,i,1) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) b_mat(ipoint,k,i,1) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1))
b_mat(ipoint,k,i,2) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) b_mat(ipoint,k,i,2) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2))
b_mat(ipoint,k,i,3) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)) b_mat(ipoint,k,i,3) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3))
@ -173,19 +152,19 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
ac_mat = 0.d0 ac_mat = 0.d0
do m = 1, 3 do m = 1, 3
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
, int2_grad1_u12_ao_test(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid & , int2_grad1_u12_ao_test(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
, 1.d0, ac_mat, ao_num*ao_num) , 1.d0, ac_mat, ao_num*ao_num)
enddo enddo
deallocate(b_mat) deallocate(b_mat)
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l) & !$OMP PRIVATE (i, j, k, l) &
!$OMP SHARED (ac_mat, tc_grad_and_lapl_ao_test, ao_num) !$OMP SHARED (ac_mat, tc_grad_and_lapl_ao_test, ao_num)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do j = 1, ao_num do j = 1, ao_num
@ -199,29 +178,23 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
deallocate(ac_mat) deallocate(ac_mat)
endif endif
if(write_tc_integ) then if(write_tc_integ.and.mpi_master) then
open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao_test', action="write") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_and_lapl_ao_test', action="write")
do i = 1, ao_num call ezfio_set_work_empty(.False.)
do j = 1, ao_num write(11) tc_grad_and_lapl_ao_test
do k = 1, ao_num
do l = 1, ao_num
write(11) tc_grad_and_lapl_ao_test(l,k,j,i)
enddo
enddo
enddo
enddo
close(11) close(11)
call ezfio_set_tc_keywords_io_tc_integ('Read')
endif endif
call wall_time(time1) call wall_time(time1)
print*, ' Wall time for tc_grad_and_lapl_ao_test = ', time1 - time0 print*, ' Wall time for tc_grad_and_lapl_ao_test = ', time1 - time0
END_PROVIDER END_PROVIDER
! --- ! ---

34
src/tc_keywords/EZFIO.cfg
View File

@ -6,7 +6,7 @@ default: False
[comp_left_eigv] [comp_left_eigv]
type: logical type: logical
doc: If |true|, computes also the left-eigenvector doc: If |true|, computes also the left-eigenvector
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
@ -14,7 +14,7 @@ default: False
type: logical type: logical
doc: If |true|, three-body terms are included doc: If |true|, three-body terms are included
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: True default: True
[pure_three_body_h_tc] [pure_three_body_h_tc]
type: logical type: logical
@ -30,13 +30,13 @@ default: False
[core_tc_op] [core_tc_op]
type: logical type: logical
doc: If |true|, takes the usual Hamiltonian for core orbitals (assumed to be doubly occupied) doc: If |true|, takes the usual Hamiltonian for core orbitals (assumed to be doubly occupied)
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
[full_tc_h_solver] [full_tc_h_solver]
type: logical type: logical
doc: If |true|, you diagonalize the full TC H matrix doc: If |true|, you diagonalize the full TC H matrix
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
@ -60,11 +60,11 @@ default: 0.000005
[thresh_psi_r_norm] [thresh_psi_r_norm]
type: logical type: logical
doc: If |true|, you prune the WF to compute the PT1 coef based on the norm. If False, the pruning is done through the amplitude on the right-coefficient. doc: If |true|, you prune the WF to compute the PT1 coef based on the norm. If False, the pruning is done through the amplitude on the right-coefficient.
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
[state_following_tc] [state_following_tc]
type: logical type: logical
doc: If |true|, the states are re-ordered to match the input states doc: If |true|, the states are re-ordered to match the input states
default: False default: False
@ -78,7 +78,7 @@ default: True
[symetric_fock_tc] [symetric_fock_tc]
type: logical type: logical
doc: If |true|, using F+F^t as Fock TC doc: If |true|, using F+F^t as Fock TC
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
@ -126,7 +126,7 @@ default: 1.e-6
[maxovl_tc] [maxovl_tc]
type: logical type: logical
doc: If |true|, maximize the overlap between orthogonalized left- and right eigenvectors doc: If |true|, maximize the overlap between orthogonalized left- and right eigenvectors
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
@ -152,7 +152,7 @@ default: 0.
type: character*(32) type: character*(32)
doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS] doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: Simple default: DIIS
[im_thresh_tcscf] [im_thresh_tcscf]
type: Threshold type: Threshold
@ -180,21 +180,15 @@ default: 1.e-6
[var_tc] [var_tc]
type: logical type: logical
doc: If |true|, use VAR-TC doc: If |true|, use VAR-TC
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
[read_tc_integ] [io_tc_integ]
type: logical type: Disk_access
doc: If |true|, read integrals: int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao doc: Read/Write integrals int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: None
[write_tc_integ]
type: logical
doc: If |true|, write integrals: int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao
interface: ezfio,provider,ocaml
default: False
[debug_tc_pt2] [debug_tc_pt2]
type: integer type: integer