PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 20:34:58 +01:00
Code Releases Activity

Vectorizing integrals

Browse Source
This commit is contained in:
Anthony Scemama 2022-11-20 22:29:58 +01:00
parent 76eeade4d5
commit 768e1ac5d8
3 changed files with 104 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f
View File

@ -90,7 +90,7 @@ subroutine overlap_gauss_r12_v(D_center_,delta,A_center,B_center,power_A,power_B
thr = 1.d-10
d(:) = 0
maxab = maxval(d(1:3))
maxab = maxval(power_A(1:3))
double precision, allocatable :: D_center(:,:)
allocate(D_center(3,n_points))

135
src/utils/integration.irp.f
View File

@ -161,52 +161,82 @@ subroutine give_explicit_poly_and_gaussian_v(P_new, ldp, P_center,p,fact_k,iorde
call gaussian_product_v(alpha,A_center,beta,B_center,fact_k,p,P_center,n_points)
!TODO TRANSP
double precision, allocatable :: P_a_(:,:,:), P_b_(:,:,:), A_center_(:,:), P_center_(:,:), P_new_(:,:,:)
allocate(A_center_(n_points,3), P_center_(n_points,3), P_new_(n_points,0:ldp,3))
A_center_(1:n_points,1) = A_center(1,1:n_points)
A_center_(1:n_points,2) = A_center(2,1:n_points)
A_center_(1:n_points,3) = A_center(3,1:n_points)
P_center_(1:n_points,1) = P_center(1,1:n_points)
P_center_(1:n_points,2) = P_center(2,1:n_points)
P_center_(1:n_points,3) = P_center(3,1:n_points)
if ( ior(ior(b(1),b(2)),b(3)) == 0 ) then ! b == (0,0,0)
lda = maxval(a)
ldb = 0
allocate(P_a(0:lda,3,n_points),P_b(0:0,3,n_points))
allocate(P_a_(n_points,0:lda,3), P_b_(n_points,0:0,3))
call recentered_poly2_v0(P_a,lda,A_center,P_center,a,P_b,B_center,P_center,n_points)
call recentered_poly2_v0(P_a_,lda,A_center_,P_center_,a,P_b_,B_center,P_center_,n_points)
iorder(1:3) = a(1:3)
do ipoint=1,n_points
do xyz=1,3
P_new(0,xyz,ipoint) = P_b(0,xyz,ipoint) * P_a(0,xyz,ipoint)
P_new_(ipoint,0,xyz) = P_a_(ipoint,0,xyz) * P_b_(ipoint,0,xyz)
do i=1,a(xyz)
P_new(i,xyz,ipoint) = P_new(i,xyz,ipoint) + P_b(0,xyz,ipoint) * P_a(i,xyz,ipoint)
P_new_(ipoint,i,xyz) = P_new_(ipoint,i,xyz) + P_b_(ipoint,0,xyz) * P_a_(ipoint,i,xyz)
enddo
enddo
enddo
do ipoint=1,n_points
do i=0,ldp
P_new(i,1,ipoint) = P_new_(ipoint,i,1)
P_new(i,2,ipoint) = P_new_(ipoint,i,2)
P_new(i,3,ipoint) = P_new_(ipoint,i,3)
enddo
enddo
return
endif
lda = maxval(a)
ldb = maxval(b)
allocate(P_a(0:lda,3,n_points),P_b(0:ldb,3,n_points))
call recentered_poly2_v(P_a,lda,A_center,P_center,a,P_b,ldb,B_center,P_center,b,n_points)
allocate(P_a_(n_points,0:lda,3), P_b_(n_points,0:ldb,3))
call recentered_poly2_v(P_a_,lda,A_center_,P_center_,a,P_b_,ldb,B_center,P_center_,b,n_points)
iorder(1:3) = a(1:3) + b(1:3)
do xyz=1,3
if (b(xyz) == 0) then
do ipoint=1,n_points
P_new(0,xyz,ipoint) = P_b(0,xyz,ipoint) * P_a(0,xyz,ipoint)
P_new_(ipoint,0,xyz) = P_a_(ipoint,0,xyz) * P_b_(ipoint,0,xyz)
do i=1,a(xyz)
P_new(i,xyz,ipoint) = P_new(i,xyz,ipoint) + P_b(0,xyz,ipoint) * P_a(i,xyz,ipoint)
P_new_(ipoint,i,xyz) = P_new_(ipoint,i,xyz) + P_b_(ipoint,0,xyz) * P_a_(ipoint,i,xyz)
enddo
enddo
else
do ipoint=1,n_points
do i=0,iorder(xyz)
P_new(i,xyz,ipoint) = 0.d0
do i=0,iorder(xyz)
do ipoint=1,n_points
P_new_(ipoint,i,xyz) = 0.d0
enddo
n_new=0
call multiply_poly(P_a(0,xyz,ipoint),a(xyz),P_b(0,xyz,ipoint),b(xyz),P_new(0,xyz,ipoint),n_new)
enddo
call multiply_poly_v(P_a_(1,0,xyz), a(xyz),P_b_(1,0,xyz),b(xyz),P_new_(1,0,xyz),ldp,n_points)
endif
enddo
do ipoint=1,n_points
do i=0,ldp
P_new(i,1,ipoint) = P_new_(ipoint,i,1)
P_new(i,2,ipoint) = P_new_(ipoint,i,2)
P_new(i,3,ipoint) = P_new_(ipoint,i,3)
enddo
enddo
end
!-
@ -487,6 +517,45 @@ subroutine multiply_poly(b,nb,c,nc,d,nd)
end
subroutine multiply_poly_v(b,nb,c,nc,d,nd,n_points)
implicit none
BEGIN_DOC
! Multiply pairs of polynomials
! D(t) += B(t)*C(t)
END_DOC
integer, intent(in) :: nb, nc, n_points
integer, intent(in) :: nd
double precision, intent(in) :: b(n_points,0:nb), c(n_points,0:nc)
double precision, intent(inout) :: d(n_points,0:nd)
integer :: ib, ic, id, k, ipoint
if (nd < nb+nc) then
print *, nd, nb, nc
print *, irp_here, ': nd < nb+nc'
stop 1
endif
do ic = 0,nc
do ipoint=1, n_points
d(ipoint,ic) = d(ipoint,ic) + c(ipoint,ic) * b(ipoint,0)
enddo
enddo
do ib=1,nb
do ipoint=1, n_points
d(ipoint, ib) = d(ipoint, ib) + c(ipoint,0) * b(ipoint, ib)
enddo
do ic = 1,nc
do ipoint=1, n_points
d(ipoint, ib+ic) = d(ipoint, ib+ic) + c(ipoint,ic) * b(ipoint, ib)
enddo
enddo
enddo
end
subroutine add_poly(b,nb,c,nc,d,nd)
implicit none
BEGIN_DOC
@ -593,8 +662,8 @@ subroutine recentered_poly2_v(P_new,lda,x_A,x_P,a,P_new2,ldb,x_B,x_Q,b,n_points)
! Recenter two polynomials
END_DOC
integer, intent(in) :: a(3),b(3), n_points, lda, ldb
double precision, intent(in) :: x_A(3,n_points),x_P(3,n_points),x_B(3),x_Q(3,n_points)
double precision, intent(out) :: P_new(0:lda,3,n_points),P_new2(0:ldb,3,n_points)
double precision, intent(in) :: x_A(n_points,3),x_P(n_points,3),x_B(3),x_Q(n_points,3)
double precision, intent(out) :: P_new(n_points,0:lda,3),P_new2(n_points,0:ldb,3)
double precision :: binom_func
integer :: i,j,k,l, minab(3), maxab(3),ipoint, xyz
double precision, allocatable :: pows_a(:,:), pows_b(:,:)
@ -610,9 +679,9 @@ subroutine recentered_poly2_v(P_new,lda,x_A,x_P,a,P_new2,ldb,x_B,x_Q,b,n_points)
if ((a(xyz)<0).or.(b(xyz)<0) ) cycle
do ipoint=1,n_points
pows_a(ipoint,0) = 1.d0
pows_a(ipoint,1) = (x_P(xyz,ipoint) - x_A(xyz,ipoint))
pows_a(ipoint,1) = (x_P(ipoint,xyz) - x_A(ipoint,xyz))
pows_b(ipoint,0) = 1.d0
pows_b(ipoint,1) = (x_Q(xyz,ipoint) - x_B(xyz))
pows_b(ipoint,1) = (x_Q(ipoint,xyz) - x_B(xyz))
enddo
do i = 2,maxab(xyz)
do ipoint=1,n_points
@ -621,39 +690,39 @@ subroutine recentered_poly2_v(P_new,lda,x_A,x_P,a,P_new2,ldb,x_B,x_Q,b,n_points)
enddo
enddo
do ipoint=1,n_points
P_new (0,xyz,ipoint) = pows_a(ipoint,a(xyz))
P_new2(0,xyz,ipoint) = pows_b(ipoint,b(xyz))
P_new (ipoint,0,xyz) = pows_a(ipoint,a(xyz))
P_new2(ipoint,0,xyz) = pows_b(ipoint,b(xyz))
enddo
do i = 1,min(minab(xyz),20)
fa = binom_transp(a(xyz)-i,a(xyz))
fb = binom_transp(b(xyz)-i,b(xyz))
do ipoint=1,n_points
P_new (i,xyz,ipoint) = fa * pows_a(ipoint,a(xyz)-i)
P_new2(i,xyz,ipoint) = fb * pows_b(ipoint,b(xyz)-i)
P_new (ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i)
P_new2(ipoint,i,xyz) = fb * pows_b(ipoint,b(xyz)-i)
enddo
enddo
do i = minab(xyz)+1,min(a(xyz),20)
fa = binom_transp(a(xyz)-i,a(xyz))
do ipoint=1,n_points
P_new (i,xyz,ipoint) = fa * pows_a(ipoint,a(xyz)-i)
P_new (ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i)
enddo
enddo
do i = minab(xyz)+1,min(b(xyz),20)
fb = binom_transp(b(xyz)-i,b(xyz))
do ipoint=1,n_points
P_new2(i,xyz,ipoint) = fb * pows_b(ipoint,b(xyz)-i)
P_new2(ipoint,i,xyz) = fb * pows_b(ipoint,b(xyz)-i)
enddo
enddo
do i = 21,a(xyz)
fa = binom_func(a(xyz),a(xyz)-i)
do ipoint=1,n_points
P_new (i,xyz,ipoint) = fa * pows_a(ipoint,a(xyz)-i)
P_new (ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i)
enddo
enddo
do i = 21,b(xyz)
fb = binom_func(b(xyz),b(xyz)-i)
do ipoint=1,n_points
P_new2(i,xyz,ipoint) = fb * pows_b(ipoint,b(xyz)-i)
P_new2(ipoint,i,xyz) = fb * pows_b(ipoint,b(xyz)-i)
enddo
enddo
enddo
@ -666,8 +735,8 @@ subroutine recentered_poly2_v0(P_new,lda,x_A,x_P,a,P_new2,x_B,x_Q,n_points)
! Recenter two polynomials. Special case for b=(0,0,0)
END_DOC
integer, intent(in) :: a(3), n_points, lda
double precision, intent(in) :: x_A(3,n_points),x_P(3,n_points),x_B(3),x_Q(3,n_points)
double precision, intent(out) :: P_new(0:lda,3,n_points),P_new2(3,n_points)
double precision, intent(in) :: x_A(n_points,3),x_P(n_points,3),x_B(3),x_Q(n_points,3)
double precision, intent(out) :: P_new(n_points,0:lda,3),P_new2(n_points,3)
double precision :: binom_func
integer :: i,j,k,l, xyz, ipoint, maxab(3)
double precision, allocatable :: pows_a(:,:), pows_b(:,:)
@ -681,9 +750,9 @@ subroutine recentered_poly2_v0(P_new,lda,x_A,x_P,a,P_new2,x_B,x_Q,n_points)
if (a(xyz)<0) cycle
do ipoint=1,n_points
pows_a(ipoint,0) = 1.d0
pows_a(ipoint,1) = (x_P(xyz,ipoint) - x_A(xyz,ipoint))
pows_a(ipoint,1) = (x_P(ipoint,xyz) - x_A(ipoint,xyz))
pows_b(ipoint,0) = 1.d0
pows_b(ipoint,1) = (x_Q(xyz,ipoint) - x_B(xyz))
pows_b(ipoint,1) = (x_Q(ipoint,xyz) - x_B(xyz))
enddo
do i = 2,maxab(xyz)
do ipoint=1,n_points
@ -692,19 +761,19 @@ subroutine recentered_poly2_v0(P_new,lda,x_A,x_P,a,P_new2,x_B,x_Q,n_points)
enddo
enddo
do ipoint=1,n_points
P_new (0,xyz,ipoint) = pows_a(ipoint,a(xyz))
P_new2(xyz,ipoint) = pows_b(ipoint,0)
P_new (ipoint,0,xyz) = pows_a(ipoint,a(xyz))
P_new2(ipoint,xyz) = pows_b(ipoint,0)
enddo
do i = 1,min(a(xyz),20)
fa = binom_transp(a(xyz)-i,a(xyz))
do ipoint=1,n_points
P_new (i,xyz,ipoint) = fa * pows_a(ipoint,a(xyz)-i)
P_new (ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i)
enddo
enddo
do i = 21,a(xyz)
fa = binom_func(a(xyz),a(xyz)-i)
do ipoint=1,n_points
P_new (i,xyz,ipoint) = fa * pows_a(ipoint,a(xyz)-i)
P_new (ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i)
enddo
enddo

2
src/utils/one_e_integration.irp.f
View File

@ -176,7 +176,7 @@ subroutine overlap_gaussian_xyz_v(A_center,B_center,alpha,beta,power_A,&
integer :: nmax
double precision :: F_integral
ldp = max_dim
ldp = maxval( power_A(1:3) + power_B(1:3) )
allocate(P_new(0:ldp,3,n_points), P_center(3,n_points), fact_p(n_points), &
fact_pp(n_points), pp(n_points))