From 6971bf186cf020ce66d0bac091d06ae850bd803f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 1 Jun 2023 21:42:02 +0200 Subject: [PATCH] Accelerated multiply_poly --- src/ao_one_e_ints/pot_ao_ints.irp.f | 21 +- src/ao_two_e_ints/two_e_integrals.irp.f | 136 +------ src/bi_ort_ints/three_body_ints_bi_ort.irp.f | 52 +-- src/utils/integration.irp.f | 366 +++++++++++++------ 4 files changed, 305 insertions(+), 270 deletions(-) diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index 928053ad..446bf730 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -455,10 +455,12 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in) do ix=0,nx X(ix) *= dble(c) enddo - call multiply_poly(X,nx,R2x,2,d,nd) +! call multiply_poly(X,nx,R2x,2,d,nd) + call multiply_poly_c2(X,nx,R2x,d,nd) ny=0 call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,Y,ny,n_pt_in) - call multiply_poly(Y,ny,R1x,2,d,nd) +! call multiply_poly(Y,ny,R1x,2,d,nd) + call multiply_poly_c2(Y,ny,R1x,d,nd) else do ix=0,n_pt_in X(ix) = 0.d0 @@ -469,7 +471,8 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in) do ix=0,nx X(ix) *= dble(a-1) enddo - call multiply_poly(X,nx,R2x,2,d,nd) +! call multiply_poly(X,nx,R2x,2,d,nd) + call multiply_poly_c2(X,nx,R2x,d,nd) nx = nd do ix=0,n_pt_in @@ -479,10 +482,12 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in) do ix=0,nx X(ix) *= dble(c) enddo - call multiply_poly(X,nx,R2x,2,d,nd) +! call multiply_poly(X,nx,R2x,2,d,nd) + call multiply_poly_c2(X,nx,R2x,d,nd) ny=0 call I_x1_pol_mult_one_e(a-1,c,R1x,R1xp,R2x,Y,ny,n_pt_in) - call multiply_poly(Y,ny,R1x,2,d,nd) +! call multiply_poly(Y,ny,R1x,2,d,nd) + call multiply_poly_c2(Y,ny,R1x,d,nd) endif end @@ -519,7 +524,8 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim) do ix=0,nx X(ix) *= dble(c-1) enddo - call multiply_poly(X,nx,R2x,2,d,nd) +! call multiply_poly(X,nx,R2x,2,d,nd) + call multiply_poly_c2(X,nx,R2x,d,nd) ny = 0 do ix=0,dim Y(ix) = 0.d0 @@ -527,7 +533,8 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim) call I_x1_pol_mult_one_e(0,c-1,R1x,R1xp,R2x,Y,ny,dim) if(ny.ge.0)then - call multiply_poly(Y,ny,R1xp,2,d,nd) +! call multiply_poly(Y,ny,R1xp,2,d,nd) + call multiply_poly_c2(Y,ny,R1xp,d,nd) endif endif end diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index 835dc89a..85ff5bcf 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -975,18 +975,7 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt ! !DIR$ FORCEINLINE ! call multiply_poly(X,nx,B_10,2,d,nd) - if (nx >= 0) then - integer :: ib - do ib=0,nx - d(ib ) = d(ib ) + B_10(0) * X(ib) - d(ib+1) = d(ib+1) + B_10(1) * X(ib) - d(ib+2) = d(ib+2) + B_10(2) * X(ib) - enddo - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(X,nx,B_10,d,nd) nx = nd !DIR$ LOOP COUNT(8) @@ -1009,17 +998,7 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt endif ! !DIR$ FORCEINLINE ! call multiply_poly(X,nx,B_00,2,d,nd) - if (nx >= 0) then - do ib=0,nx - d(ib ) = d(ib ) + B_00(0) * X(ib) - d(ib+1) = d(ib+1) + B_00(1) * X(ib) - d(ib+2) = d(ib+2) + B_00(2) * X(ib) - enddo - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(X,nx,B_00,d,nd) endif ny=0 @@ -1038,17 +1017,7 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt ! !DIR$ FORCEINLINE ! call multiply_poly(Y,ny,C_00,2,d,nd) - if (ny >= 0) then - do ib=0,ny - d(ib ) = d(ib ) + C_00(0) * Y(ib) - d(ib+1) = d(ib+1) + C_00(1) * Y(ib) - d(ib+2) = d(ib+2) + C_00(2) * Y(ib) - enddo - - do nd = ny+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(Y,ny,C_00,d,nd) end recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) @@ -1088,18 +1057,7 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) ! !DIR$ FORCEINLINE ! call multiply_poly(X,nx,B_00,2,d,nd) - if (nx >= 0) then - integer :: ib - do ib=0,nx - d(ib ) = d(ib ) + B_00(0) * X(ib) - d(ib+1) = d(ib+1) + B_00(1) * X(ib) - d(ib+2) = d(ib+2) + B_00(2) * X(ib) - enddo - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(X,nx,B_00,d,nd) ny=0 @@ -1111,17 +1069,7 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) ! !DIR$ FORCEINLINE ! call multiply_poly(Y,ny,C_00,2,d,nd) - if (ny >= 0) then - do ib=0,ny - d(ib ) = d(ib ) + C_00(0) * Y(ib) - d(ib+1) = d(ib+1) + C_00(1) * Y(ib) - d(ib+2) = d(ib+2) + C_00(2) * Y(ib) - enddo - - do nd = ny+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(Y,ny,C_00,d,nd) end @@ -1150,18 +1098,7 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) ! !DIR$ FORCEINLINE ! call multiply_poly(X,nx,B_10,2,d,nd) - if (nx >= 0) then - integer :: ib - do ib=0,nx - d(ib ) = d(ib ) + B_10(0) * X(ib) - d(ib+1) = d(ib+1) + B_10(1) * X(ib) - d(ib+2) = d(ib+2) + B_10(2) * X(ib) - enddo - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(X,nx,B_10,d,nd) nx = nd !DIR$ LOOP COUNT(8) @@ -1181,17 +1118,7 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) ! !DIR$ FORCEINLINE ! call multiply_poly(X,nx,B_00,2,d,nd) - if (nx >= 0) then - do ib=0,nx - d(ib ) = d(ib ) + B_00(0) * X(ib) - d(ib+1) = d(ib+1) + B_00(1) * X(ib) - d(ib+2) = d(ib+2) + B_00(2) * X(ib) - enddo - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(X,nx,B_00,d,nd) ny=0 !DIR$ LOOP COUNT(8) @@ -1203,17 +1130,7 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) ! !DIR$ FORCEINLINE ! call multiply_poly(Y,ny,C_00,2,d,nd) - if (ny >= 0) then - do ib=0,ny - d(ib ) = d(ib ) + C_00(0) * Y(ib) - d(ib+1) = d(ib+1) + C_00(1) * Y(ib) - d(ib+2) = d(ib+2) + C_00(2) * Y(ib) - enddo - - do nd = ny+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(Y,ny,C_00,d,nd) end recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim) @@ -1262,18 +1179,7 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim) ! !DIR$ FORCEINLINE ! call multiply_poly(Y,ny,D_00,2,d,nd) - if (ny >= 0) then - integer :: ib - do ib=0,ny - d(ib ) = d(ib ) + D_00(0) * Y(ib) - d(ib+1) = d(ib+1) + D_00(1) * Y(ib) - d(ib+2) = d(ib+2) + D_00(2) * Y(ib) - enddo - - do nd = ny+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(Y,ny,D_00,d,nd) return @@ -1293,17 +1199,7 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim) ! !DIR$ FORCEINLINE ! call multiply_poly(X,nx,B_01,2,d,nd) - if (nx >= 0) then - do ib=0,nx - d(ib ) = d(ib ) + B_01(0) * X(ib) - d(ib+1) = d(ib+1) + B_01(1) * X(ib) - d(ib+2) = d(ib+2) + B_01(2) * X(ib) - enddo - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(X,nx,B_01,d,nd) ny = 0 !DIR$ LOOP COUNT(6) @@ -1314,17 +1210,7 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim) ! !DIR$ FORCEINLINE ! call multiply_poly(Y,ny,D_00,2,d,nd) - if (ny >= 0) then - do ib=0,ny - d(ib ) = d(ib ) + D_00(0) * Y(ib) - d(ib+1) = d(ib+1) + D_00(1) * Y(ib) - d(ib+2) = d(ib+2) + D_00(2) * Y(ib) - enddo - - do nd = ny+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - endif + call multiply_poly_c2(Y,ny,D_00,d,nd) end select end diff --git a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f index e8b56307..a72cd682 100644 --- a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f +++ b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f @@ -4,7 +4,7 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num, mo_num)] BEGIN_DOC -! matrix element of the -L three-body operator +! matrix element of the -L three-body operator ! ! notice the -1 sign: in this way three_body_ints_bi_ort can be directly used to compute Slater rules :) END_DOC @@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n implicit none integer :: i, j, k, l, m, n double precision :: integral, wall1, wall0 - character*(128) :: name_file + character*(128) :: name_file three_body_ints_bi_ort = 0.d0 print *, ' Providing the three_body_ints_bi_ort ...' @@ -27,12 +27,12 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n ! call read_array_6_index_tensor(mo_num,three_body_ints_bi_ort,name_file) ! else - !provide x_W_ki_bi_ortho_erf_rk + !provide x_W_ki_bi_ortho_erf_rk provide mos_r_in_r_array_transp mos_l_in_r_array_transp !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,n,integral) & + !$OMP PRIVATE (i,j,k,l,m,n,integral) & !$OMP SHARED (mo_num,three_body_ints_bi_ort) !$OMP DO SCHEDULE (dynamic) do i = 1, mo_num @@ -43,7 +43,7 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n do n = 1, mo_num call give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) - three_body_ints_bi_ort(n,l,k,m,j,i) = -1.d0 * integral + three_body_ints_bi_ort(n,l,k,m,j,i) = -1.d0 * integral enddo enddo enddo @@ -63,7 +63,7 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n ! call ezfio_set_three_body_ints_bi_ort_io_three_body_ints_bi_ort("Read") ! endif -END_PROVIDER +END_PROVIDER ! --- @@ -71,7 +71,7 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) BEGIN_DOC ! - ! < n l k | -L | m j i > with a BI-ORTHONORMAL MOLECULAR ORBITALS + ! < n l k | -L | m j i > with a BI-ORTHONORMAL MOLECULAR ORBITALS ! END_DOC @@ -79,28 +79,30 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) integer, intent(in) :: n, l, k, m, j, i double precision, intent(out) :: integral integer :: ipoint - double precision :: weight + double precision :: weight, tmp PROVIDE mo_l_coef mo_r_coef PROVIDE int2_grad1_u12_bimo_t integral = 0.d0 do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & + tmp = mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & * ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,l,j) & + int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,l,j) & + int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,l,j) ) - integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & + + tmp = tmp + mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & * ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,k,i) & + int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,k,i) & + int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,k,i) ) - integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & + + tmp = tmp + mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & * ( int2_grad1_u12_bimo_t(ipoint,1,l,j) * int2_grad1_u12_bimo_t(ipoint,1,k,i) & + int2_grad1_u12_bimo_t(ipoint,2,l,j) * int2_grad1_u12_bimo_t(ipoint,2,k,i) & + int2_grad1_u12_bimo_t(ipoint,3,l,j) * int2_grad1_u12_bimo_t(ipoint,3,k,i) ) + integral = integral + tmp * final_weight_at_r_vector(ipoint) enddo end subroutine give_integrals_3_body_bi_ort @@ -111,7 +113,7 @@ subroutine give_integrals_3_body_bi_ort_old(n, l, k, m, j, i, integral) BEGIN_DOC ! - ! < n l k | -L | m j i > with a BI-ORTHONORMAL MOLECULAR ORBITALS + ! < n l k | -L | m j i > with a BI-ORTHONORMAL MOLECULAR ORBITALS ! END_DOC @@ -123,13 +125,13 @@ subroutine give_integrals_3_body_bi_ort_old(n, l, k, m, j, i, integral) integral = 0.d0 do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) + weight = final_weight_at_r_vector(ipoint) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & +! integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & ! * ( x_W_ki_bi_ortho_erf_rk(ipoint,1,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,1,l,j) & ! + x_W_ki_bi_ortho_erf_rk(ipoint,2,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,2,l,j) & ! + x_W_ki_bi_ortho_erf_rk(ipoint,3,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,3,l,j) ) -! integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & +! integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & ! * ( x_W_ki_bi_ortho_erf_rk(ipoint,1,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,1,k,i) & ! + x_W_ki_bi_ortho_erf_rk(ipoint,2,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,2,k,i) & ! + x_W_ki_bi_ortho_erf_rk(ipoint,3,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,3,k,i) ) @@ -138,11 +140,11 @@ subroutine give_integrals_3_body_bi_ort_old(n, l, k, m, j, i, integral) ! + x_W_ki_bi_ortho_erf_rk(ipoint,2,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,2,k,i) & ! + x_W_ki_bi_ortho_erf_rk(ipoint,3,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,3,k,i) ) -! integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & +! integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & ! * ( int2_grad1_u12_bimo(1,n,m,ipoint) * int2_grad1_u12_bimo(1,l,j,ipoint) & ! + int2_grad1_u12_bimo(2,n,m,ipoint) * int2_grad1_u12_bimo(2,l,j,ipoint) & ! + int2_grad1_u12_bimo(3,n,m,ipoint) * int2_grad1_u12_bimo(3,l,j,ipoint) ) -! integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & +! integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & ! * ( int2_grad1_u12_bimo(1,n,m,ipoint) * int2_grad1_u12_bimo(1,k,i,ipoint) & ! + int2_grad1_u12_bimo(2,n,m,ipoint) * int2_grad1_u12_bimo(2,k,i,ipoint) & ! + int2_grad1_u12_bimo(3,n,m,ipoint) * int2_grad1_u12_bimo(3,k,i,ipoint) ) @@ -151,13 +153,13 @@ subroutine give_integrals_3_body_bi_ort_old(n, l, k, m, j, i, integral) ! + int2_grad1_u12_bimo(2,l,j,ipoint) * int2_grad1_u12_bimo(2,k,i,ipoint) & ! + int2_grad1_u12_bimo(3,l,j,ipoint) * int2_grad1_u12_bimo(3,k,i,ipoint) ) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & + integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & * ( int2_grad1_u12_bimo_transp(n,m,1,ipoint) * int2_grad1_u12_bimo_transp(l,j,1,ipoint) & + int2_grad1_u12_bimo_transp(n,m,2,ipoint) * int2_grad1_u12_bimo_transp(l,j,2,ipoint) & + int2_grad1_u12_bimo_transp(n,m,3,ipoint) * int2_grad1_u12_bimo_transp(l,j,3,ipoint) ) - integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & + integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & * ( int2_grad1_u12_bimo_transp(n,m,1,ipoint) * int2_grad1_u12_bimo_transp(k,i,1,ipoint) & + int2_grad1_u12_bimo_transp(n,m,2,ipoint) * int2_grad1_u12_bimo_transp(k,i,2,ipoint) & + int2_grad1_u12_bimo_transp(n,m,3,ipoint) * int2_grad1_u12_bimo_transp(k,i,3,ipoint) ) @@ -176,7 +178,7 @@ subroutine give_integrals_3_body_bi_ort_ao(n, l, k, m, j, i, integral) BEGIN_DOC ! - ! < n l k | -L | m j i > with a BI-ORTHONORMAL ATOMIC ORBITALS + ! < n l k | -L | m j i > with a BI-ORTHONORMAL ATOMIC ORBITALS ! END_DOC @@ -188,13 +190,13 @@ subroutine give_integrals_3_body_bi_ort_ao(n, l, k, m, j, i, integral) integral = 0.d0 do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) + weight = final_weight_at_r_vector(ipoint) - integral += weight * aos_in_r_array_transp(ipoint,k) * aos_in_r_array_transp(ipoint,i) & + integral += weight * aos_in_r_array_transp(ipoint,k) * aos_in_r_array_transp(ipoint,i) & * ( int2_grad1_u12_ao_t(ipoint,1,n,m) * int2_grad1_u12_ao_t(ipoint,1,l,j) & + int2_grad1_u12_ao_t(ipoint,2,n,m) * int2_grad1_u12_ao_t(ipoint,2,l,j) & + int2_grad1_u12_ao_t(ipoint,3,n,m) * int2_grad1_u12_ao_t(ipoint,3,l,j) ) - integral += weight * aos_in_r_array_transp(ipoint,l) * aos_in_r_array_transp(ipoint,j) & + integral += weight * aos_in_r_array_transp(ipoint,l) * aos_in_r_array_transp(ipoint,j) & * ( int2_grad1_u12_ao_t(ipoint,1,n,m) * int2_grad1_u12_ao_t(ipoint,1,k,i) & + int2_grad1_u12_ao_t(ipoint,2,n,m) * int2_grad1_u12_ao_t(ipoint,2,k,i) & + int2_grad1_u12_ao_t(ipoint,3,n,m) * int2_grad1_u12_ao_t(ipoint,3,k,i) ) diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f index b60e3bc1..21179dac 100644 --- a/src/utils/integration.irp.f +++ b/src/utils/integration.irp.f @@ -56,7 +56,7 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha, ! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 ) ! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 ) ! - ! WARNING ::: IF fact_k is too smal then: + ! WARNING ::: IF fact_k is too smal then: ! returns a "s" function centered in zero ! with an inifinite exponent and a zero polynom coef END_DOC @@ -86,7 +86,7 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha, !DIR$ FORCEINLINE call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center) if (fact_k < thresh) then - ! IF fact_k is too smal then: + ! IF fact_k is too smal then: ! returns a "s" function centered in zero ! with an inifinite exponent and a zero polynom coef P_center = 0.d0 @@ -468,114 +468,6 @@ end subroutine -subroutine multiply_poly_0c(b,c,nc,d,nd) - implicit none - BEGIN_DOC - ! Multiply two polynomials - ! D(t) += B(t)*C(t) - END_DOC - - integer, intent(in) :: nc - integer, intent(out) :: nd - double precision, intent(in) :: b(0:0), c(0:nc) - double precision, intent(inout) :: d(0:0+nc) - - integer :: ic - - do ic = 0,nc - d(ic) = d(ic) + c(ic) * b(0) - enddo - - do nd = nc,0,-1 - if (d(nd) /= 0.d0) exit - enddo - -end - -subroutine multiply_poly_1c(b,c,nc,d,nd) - implicit none - BEGIN_DOC - ! Multiply two polynomials - ! D(t) += B(t)*C(t) - END_DOC - - integer, intent(in) :: nc - integer, intent(out) :: nd - double precision, intent(in) :: b(0:1), c(0:nc) - double precision, intent(inout) :: d(0:1+nc) - - integer :: ic, id - if(nc < 0) return - - do ic = 0,nc - d( ic) = d( ic) + c(ic) * b(0) - d(1+ic) = d(1+ic) + c(ic) * b(1) - enddo - - do nd = nc+1,0,-1 - if (d(nd) /= 0.d0) exit - enddo - -end - - -subroutine multiply_poly_2c(b,c,nc,d,nd) - implicit none - BEGIN_DOC - ! Multiply two polynomials - ! D(t) += B(t)*C(t) - END_DOC - - integer, intent(in) :: nc - integer, intent(out) :: nd - double precision, intent(in) :: b(0:2), c(0:nc) - double precision, intent(inout) :: d(0:2+nc) - - integer :: ic, id, k - if (nc <0) return - - do ic = 0,nc - d( ic) = d( ic) + c(ic) * b(0) - d(1+ic) = d(1+ic) + c(ic) * b(1) - d(2+ic) = d(2+ic) + c(ic) * b(2) - enddo - - do nd = nc+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - -end - -subroutine multiply_poly_3c(b,c,nc,d,nd) - implicit none - BEGIN_DOC - ! Multiply two polynomials - ! D(t) += B(t)*C(t) - END_DOC - - integer, intent(in) :: nc - integer, intent(out) :: nd - double precision, intent(in) :: b(0:3), c(0:nc) - double precision, intent(inout) :: d(0:3+nc) - - integer :: ic, id - if (nc <0) return - - do ic = 1,nc - d( ic) = d(1+ic) + c(ic) * b(0) - d(1+ic) = d(1+ic) + c(ic) * b(1) - d(2+ic) = d(1+ic) + c(ic) * b(2) - d(3+ic) = d(1+ic) + c(ic) * b(3) - enddo - - do nd = nc+3,0,-1 - if (d(nd) /= 0.d0) exit - enddo - -end - - - subroutine multiply_poly(b,nb,c,nc,d,nd) implicit none BEGIN_DOC @@ -604,6 +496,254 @@ subroutine multiply_poly(b,nb,c,nc,d,nd) end + +subroutine multiply_poly_b0(b,c,nc,d,nd) + implicit none + BEGIN_DOC + ! Multiply two polynomials + ! D(t) += B(t)*C(t) + END_DOC + + integer, intent(in) :: nc + integer, intent(out) :: nd + double precision, intent(in) :: b(0:0), c(0:nc) + double precision, intent(inout) :: d(0:nc) + + integer :: ndtmp + integer :: ic, id, k + if(nc < 0) return !False if nc>=0 + + do ic = 0,nc + d(ic) = d(ic) + c(ic) * b(0) + enddo + + do nd = nc,0,-1 + if (d(nd) /= 0.d0) exit + enddo + +end + +subroutine multiply_poly_b1(b,c,nc,d,nd) + implicit none + BEGIN_DOC + ! Multiply two polynomials + ! D(t) += B(t)*C(t) + END_DOC + + integer, intent(in) :: nc + integer, intent(out) :: nd + double precision, intent(in) :: b(0:1), c(0:nc) + double precision, intent(inout) :: d(0:1+nc) + + integer :: ndtmp + integer :: ib, ic, id, k + if(nc < 0) return !False if nc>=0 + + + select case (nc) + case (0) + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(0) * b(1) + + case (1) + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(0) * b(1) + c(1) * b(0) + d(2) = d(2) + c(1) * b(1) + + case default + d(0) = d(0) + c(0) * b(0) + do ic = 1,nc + d(ic) = d(ic) + c(ic) * b(0) + c(ic-1) * b(1) + enddo + d(nc+1) = d(nc+1) + c(nc) * b(1) + + end select + + do nd = 1+nc,0,-1 + if (d(nd) /= 0.d0) exit + enddo + +end + + +subroutine multiply_poly_b2(b,c,nc,d,nd) + implicit none + BEGIN_DOC + ! Multiply two polynomials + ! D(t) += B(t)*C(t) + END_DOC + + integer, intent(in) :: nc + integer, intent(out) :: nd + double precision, intent(in) :: b(0:2), c(0:nc) + double precision, intent(inout) :: d(0:2+nc) + + integer :: ndtmp + integer :: ib, ic, id, k + if(nc < 0) return !False if nc>=0 + + select case (nc) + case (0) + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(0) * b(1) + d(2) = d(2) + c(0) * b(2) + + case (1) + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(0) * b(1) + c(1) * b(0) + d(2) = d(2) + c(0) * b(2) + c(1) * b(1) + d(3) = d(3) + c(1) * b(2) + + case (2) + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(0) * b(1) + c(1) * b(0) + d(2) = d(2) + c(0) * b(2) + c(1) * b(1) + c(2) * b(0) + d(3) = d(3) + c(2) * b(1) + c(1) * b(2) + d(4) = d(4) + c(2) * b(2) + + case default + + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(0) * b(1) + c(1) * b(0) + do ic = 2,nc + d(ic) = d(ic) + c(ic) * b(0) + c(ic-1) * b(1) + c(ic-2) * b(2) + enddo + d(nc+1) = d(nc+1) + c(nc) * b(1) + c(nc-1) * b(2) + d(nc+2) = d(nc+2) + c(nc) * b(2) + + end select + + do nd = 2+nc,0,-1 + if (d(nd) /= 0.d0) exit + enddo + +end + + +subroutine multiply_poly_c0(b,nb,c,d,nd) + implicit none + BEGIN_DOC + ! Multiply two polynomials + ! D(t) += B(t)*C(t) + END_DOC + + integer, intent(in) :: nb + integer, intent(out) :: nd + double precision, intent(in) :: b(0:nb), c(0:0) + double precision, intent(inout) :: d(0:nb) + + integer :: ndtmp + integer :: ib, ic, id, k + if(nb < 0) return !False if nb>=0 + + do ib=0,nb + d(ib) = d(ib) + c(0) * b(ib) + enddo + + do nd = nb,0,-1 + if (d(nd) /= 0.d0) exit + enddo + +end + + +subroutine multiply_poly_c1(b,nb,c,d,nd) + implicit none + BEGIN_DOC + ! Multiply two polynomials + ! D(t) += B(t)*C(t) + END_DOC + + integer, intent(in) :: nb + integer, intent(out) :: nd + double precision, intent(in) :: b(0:nb), c(0:1) + double precision, intent(inout) :: d(0:nb+1) + + integer :: ndtmp + integer :: ib, ic, id, k + if(nb < 0) return !False if nb>=0 + + select case (nb) + case (0) + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(1) * b(0) + + case (1) + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(0) * b(1) + c(1) * b(0) + d(2) = d(2) + c(1) * b(1) + + case default + d(0) = d(0) + c(0) * b(0) + do ib=1,nb + d(ib) = d(ib) + c(0) * b(ib) + c(1) * b(ib-1) + enddo + d(nb+1) = d(nb+1) + c(1) * b(nb) + + end select + + do nd = nb+1,0,-1 + if (d(nd) /= 0.d0) exit + enddo + +end + + +subroutine multiply_poly_c2(b,nb,c,d,nd) + implicit none + BEGIN_DOC + ! Multiply two polynomials + ! D(t) += B(t)*C(t) + END_DOC + + integer, intent(in) :: nb + integer, intent(out) :: nd + double precision, intent(in) :: b(0:nb), c(0:2) + double precision, intent(inout) :: d(0:nb+2) + + integer :: ndtmp + integer :: ib, ic, id, k + if(nb < 0) return !False if nb>=0 + + select case (nb) + case (0) + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(1) * b(0) + d(2) = d(2) + c(2) * b(0) + + case (1) + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(0) * b(1) + c(1) * b(0) + d(2) = d(2) + c(1) * b(1) + c(2) * b(0) + d(3) = d(3) + c(2) * b(1) + + case (2) + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(0) * b(1) + c(1) * b(0) + d(2) = d(2) + c(0) * b(2) + c(1) * b(1) + c(2) * b(0) + d(3) = d(3) + c(1) * b(2) + c(2) * b(1) + d(4) = d(4) + c(2) * b(2) + + case default + d(0) = d(0) + c(0) * b(0) + d(1) = d(1) + c(0) * b(1) + c(1) * b(0) + do ib=2,nb + d(ib) = d(ib) + c(0) * b(ib) + c(1) * b(ib-1) + c(2) * b(ib-2) + enddo + d(nb+1) = d(nb+1) + c(1) * b(nb) + c(2) * b(nb-1) + d(nb+2) = d(nb+2) + c(2) * b(nb) + + end select + + do nd = nb+2,0,-1 + if (d(nd) /= 0.d0) exit + enddo + +end + + + + subroutine multiply_poly_v(b,nb,c,nc,d,nd,n_points) implicit none BEGIN_DOC @@ -778,11 +918,11 @@ end subroutine recentered_poly2_v subroutine recentered_poly2_v0(P_new, lda, x_A, LD_xA, x_P, a, n_points) BEGIN_DOC - ! + ! ! Recenter two polynomials. Special case for b=(0,0,0) - ! + ! ! (x - A)^a (x - B)^0 = (x - P + P - A)^a (x - Q + Q - B)^0 - ! = (x - P + P - A)^a + ! = (x - P + P - A)^a ! END_DOC