subroutine give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_k,iorder,alpha,beta,a,b,A_center,B_center,dim) BEGIN_DOC ! Transform the product of ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) ! into ! fact_k (x-x_P)^iorder(1) (y-y_P)^iorder(2) (z-z_P)^iorder(3) exp(-p(r-P)^2) END_DOC implicit none include 'constants.include.F' integer, intent(in) :: dim integer, intent(in) :: a,b ! powers : (x-xa)**a_x = (x-A(1))**a(1) double precision, intent(in) :: alpha, beta ! exponents double precision, intent(in) :: A_center ! A center double precision, intent(in) :: B_center ! B center double precision, intent(out) :: P_center ! new center double precision, intent(out) :: p ! new exponent double precision, intent(out) :: fact_k ! constant factor double precision, intent(out) :: P_new(0:max_dim) ! polynomial integer, intent(out) :: iorder ! order of the polynomials double precision :: P_a(0:max_dim), P_b(0:max_dim) integer :: n_new,i,j double precision :: p_inv,ab,d_AB !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b ! Do the gaussian product to get the new center and the new exponent P_new = 0.d0 p = alpha+beta p_inv = 1.d0/p ab = alpha * beta d_AB = (A_center - B_center) * (A_center - B_center) P_center = (alpha * A_center + beta * B_center) * p_inv if(dabs(ab*p_inv * d_AB).lt.50.d0)then fact_k = exp(-ab*p_inv * d_AB) else fact_k = 0.d0 endif ! Recenter the polynomials P_a and P_b on x !DIR$ FORCEINLINE call recentered_poly2(P_a(0),A_center,P_center,a,P_b(0),B_center,P_center,b) n_new = 0 !DIR$ FORCEINLINE call multiply_poly(P_a(0),a,P_b(0),b,P_new(0),n_new) iorder = a + b end subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,beta,a,b,A_center,B_center,dim) BEGIN_DOC ! Transforms the product of ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) ! into ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 ) ! * [ 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: ! returns a "s" function centered in zero ! with an inifinite exponent and a zero polynom coef END_DOC implicit none include 'constants.include.F' integer, intent(in) :: dim integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) double precision, intent(in) :: alpha, beta ! exponents double precision, intent(in) :: A_center(3) ! A center double precision, intent(in) :: B_center (3) ! B center double precision, intent(out) :: P_center(3) ! new center double precision, intent(out) :: p ! new exponent double precision, intent(out) :: fact_k ! constant factor double precision, intent(out) :: P_new(0:max_dim,3)! polynomial integer, intent(out) :: iorder(3) ! i_order(i) = order of the polynomials double precision :: P_a(0:max_dim,3), P_b(0:max_dim,3) integer :: n_new,i,j !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b iorder(1) = 0 iorder(2) = 0 iorder(3) = 0 P_new(0,1) = 0.d0 P_new(0,2) = 0.d0 P_new(0,3) = 0.d0 !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: ! returns a "s" function centered in zero ! with an inifinite exponent and a zero polynom coef P_center = 0.d0 p = 1.d+15 P_new = 0.d0 iorder = 0 fact_k = 0.d0 return endif !DIR$ FORCEINLINE call recentered_poly2(P_a(0,1),A_center(1),P_center(1),a(1),P_b(0,1),B_center(1),P_center(1),b(1)) iorder(1) = a(1) + b(1) do i=0,iorder(1) P_new(i,1) = 0.d0 enddo n_new=0 !DIR$ FORCEINLINE call multiply_poly(P_a(0,1),a(1),P_b(0,1),b(1),P_new(0,1),n_new) !DIR$ FORCEINLINE call recentered_poly2(P_a(0,2),A_center(2),P_center(2),a(2),P_b(0,2),B_center(2),P_center(2),b(2)) iorder(2) = a(2) + b(2) do i=0,iorder(2) P_new(i,2) = 0.d0 enddo n_new=0 !DIR$ FORCEINLINE call multiply_poly(P_a(0,2),a(2),P_b(0,2),b(2),P_new(0,2),n_new) !DIR$ FORCEINLINE call recentered_poly2(P_a(0,3),A_center(3),P_center(3),a(3),P_b(0,3),B_center(3),P_center(3),b(3)) iorder(3) = a(3) + b(3) do i=0,iorder(3) P_new(i,3) = 0.d0 enddo n_new=0 !DIR$ FORCEINLINE call multiply_poly(P_a(0,3),a(3),P_b(0,3),b(3),P_new(0,3),n_new) end subroutine give_explicit_poly_and_gaussian_v(P_new, ldp, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, LD_A, B_center, n_points) BEGIN_DOC ! Transforms the product of ! (x-x_A)^a(1) (x-x_B)^b(1) (y-y_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) ! into ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 ) ! * [ 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: ! returns a "s" function centered in zero ! with an inifinite exponent and a zero polynom coef END_DOC include 'constants.include.F' implicit none integer, intent(in) :: n_points, ldp, LD_A integer, intent(in) :: a(3), b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) double precision, intent(in) :: alpha, beta ! exponents double precision, intent(in) :: A_center(LD_A,3) ! A center double precision, intent(in) :: B_center(3) ! B center integer, intent(out) :: iorder(3) ! i_order(i) = order of the polynomials double precision, intent(out) :: P_center(n_points,3) ! new center double precision, intent(out) :: p ! new exponent double precision, intent(out) :: fact_k(n_points) ! constant factor double precision, intent(out) :: P_new(n_points,0:ldp,3) ! polynomial integer :: n_new, i, j, ipoint, lda, ldb, xyz double precision, allocatable :: P_a(:,:,:), P_b(:,:,:) call gaussian_product_v(alpha, A_center, LD_A, beta, B_center, fact_k, p, P_center, n_points) if(ior(ior(b(1), b(2)), b(3)) == 0) then ! b == (0,0,0) iorder(1:3) = a(1:3) lda = maxval(a) allocate(P_a(n_points,0:lda,3)) !ldb = 0 !allocate(P_b(n_points,0:0,3)) !call recentered_poly2_v0(P_a, lda, A_center, LD_A, P_center, a, P_b, B_center, P_center, n_points) call recentered_poly2_v0(P_a, lda, A_center, LD_A, P_center, a, n_points) do ipoint = 1, n_points do xyz = 1, 3 !P_new(ipoint,0,xyz) = P_a(ipoint,0,xyz) * P_b(ipoint,0,xyz) P_new(ipoint,0,xyz) = P_a(ipoint,0,xyz) do i = 1, a(xyz) !P_new(ipoint,i,xyz) = P_new(ipoint,i,xyz) + P_b(ipoint,0,xyz) * P_a(ipoint,i,xyz) P_new(ipoint,i,xyz) = P_a(ipoint,i,xyz) enddo enddo enddo deallocate(P_a) !deallocate(P_b) return endif lda = maxval(a) ldb = maxval(b) allocate(P_a(n_points,0:lda,3), P_b(n_points,0:ldb,3)) call recentered_poly2_v(P_a, lda, A_center, LD_A, 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(ipoint,0,xyz) = P_a(ipoint,0,xyz) * P_b(ipoint,0,xyz) P_new(ipoint,0,xyz) = P_a(ipoint,0,xyz) do i = 1, a(xyz) !P_new(ipoint,i,xyz) = P_new(ipoint,i,xyz) + P_b(ipoint,0,xyz) * P_a(ipoint,i,xyz) P_new(ipoint,i,xyz) = P_a(ipoint,i,xyz) enddo enddo else do i = 0, iorder(xyz) do ipoint = 1, n_points P_new(ipoint,i,xyz) = 0.d0 enddo 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 end subroutine give_explicit_poly_and_gaussian_v ! --- subroutine give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_k,iorder,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,dim) BEGIN_DOC ! Transforms the product of ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) ! exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) exp(-(r-Nucl_center)^2 gama ! ! into ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 ) ! * [ 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 ) END_DOC implicit none include 'constants.include.F' integer, intent(in) :: dim integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) double precision, intent(in) :: alpha, beta, gama ! exponents double precision, intent(in) :: A_center(3) ! A center double precision, intent(in) :: B_center (3) ! B center double precision, intent(in) :: Nucl_center(3) ! B center double precision, intent(out) :: P_center(3) ! new center double precision, intent(out) :: p ! new exponent double precision, intent(out) :: fact_k ! constant factor double precision, intent(out) :: P_new(0:max_dim,3)! polynomial integer , intent(out) :: iorder(3) ! i_order(i) = order of the polynomials double precision :: P_center_tmp(3) ! new center double precision :: p_tmp ! new exponent double precision :: fact_k_tmp,fact_k_bis ! constant factor double precision :: P_new_tmp(0:max_dim,3)! polynomial integer :: i,j double precision :: binom_func ! First you transform the two primitives into a sum of primitive with the same center P_center_tmp and gaussian exponent p_tmp call give_explicit_poly_and_gaussian(P_new_tmp,P_center_tmp,p_tmp,fact_k_tmp,iorder,alpha,beta,a,b,A_center,B_center,dim) ! Then you create the new gaussian from the product of the new one per the Nuclei one call gaussian_product(p_tmp,P_center_tmp,gama,Nucl_center,fact_k_bis,p,P_center) fact_k = fact_k_bis * fact_k_tmp ! Then you build the coefficient of the new polynom do i = 0, iorder(1) P_new(i,1) = 0.d0 do j = i,iorder(1) P_new(i,1) = P_new(i,1) + P_new_tmp(j,1) * binom_func(j,j-i) * (P_center(1) - P_center_tmp(1))**(j-i) enddo enddo do i = 0, iorder(2) P_new(i,2) = 0.d0 do j = i,iorder(2) P_new(i,2) = P_new(i,2) + P_new_tmp(j,2) * binom_func(j,j-i) * (P_center(2) - P_center_tmp(2))**(j-i) enddo enddo do i = 0, iorder(3) P_new(i,3) = 0.d0 do j = i,iorder(3) P_new(i,3) = P_new(i,3) + P_new_tmp(j,3) * binom_func(j,j-i) * (P_center(3) - P_center_tmp(3))**(j-i) enddo enddo end subroutine gaussian_product(a,xa,b,xb,k,p,xp) implicit none BEGIN_DOC ! Gaussian product in 1D. ! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} END_DOC double precision, intent(in) :: a,b ! Exponents double precision, intent(in) :: xa(3),xb(3) ! Centers double precision, intent(out) :: p ! New exponent double precision, intent(out) :: xp(3) ! New center double precision, intent(out) :: k ! Constant double precision :: p_inv ASSERT (a>0.) ASSERT (b>0.) double precision :: xab(3), ab !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab p = a+b p_inv = 1.d0/(a+b) ab = a*b xab(1) = xa(1)-xb(1) xab(2) = xa(2)-xb(2) xab(3) = xa(3)-xb(3) ab = ab*p_inv k = ab*(xab(1)*xab(1)+xab(2)*xab(2)+xab(3)*xab(3)) if (k > 40.d0) then k=0.d0 return endif k = dexp(-k) xp(1) = (a*xa(1)+b*xb(1))*p_inv xp(2) = (a*xa(2)+b*xb(2))*p_inv xp(3) = (a*xa(3)+b*xb(3))*p_inv end subroutine subroutine gaussian_product_v(a, xa, LD_xa, b, xb, k, p, xp, n_points) BEGIN_DOC ! ! Gaussian product in 1D. ! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} ! ! Using multiple A centers ! END_DOC implicit none integer, intent(in) :: LD_xa, n_points double precision, intent(in) :: a, b ! Exponents double precision, intent(in) :: xa(LD_xa,3), xb(3) ! Centers double precision, intent(out) :: p ! New exponent double precision, intent(out) :: xp(n_points,3) ! New center double precision, intent(out) :: k(n_points) ! Constant integer :: ipoint double precision :: p_inv double precision :: xab(3), ab, ap, bp, bpxb(3) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab ASSERT (a>0.) ASSERT (b>0.) p = a+b p_inv = 1.d0/(a+b) ab = a*b*p_inv ap = a*p_inv bp = b*p_inv bpxb(1) = bp*xb(1) bpxb(2) = bp*xb(2) bpxb(3) = bp*xb(3) do ipoint = 1, n_points xab(1) = xa(ipoint,1)-xb(1) xab(2) = xa(ipoint,2)-xb(2) xab(3) = xa(ipoint,3)-xb(3) k(ipoint) = ab*(xab(1)*xab(1)+xab(2)*xab(2)+xab(3)*xab(3)) if (k(ipoint) > 40.d0) then k(ipoint)=0.d0 xp(ipoint,1) = 0.d0 xp(ipoint,2) = 0.d0 xp(ipoint,3) = 0.d0 else k(ipoint) = dexp(-k(ipoint)) xp(ipoint,1) = ap*xa(ipoint,1)+bpxb(1) xp(ipoint,2) = ap*xa(ipoint,2)+bpxb(2) xp(ipoint,3) = ap*xa(ipoint,3)+bpxb(3) endif enddo end subroutine gaussian_product_v ! --- subroutine gaussian_product_x(a,xa,b,xb,k,p,xp) implicit none BEGIN_DOC ! Gaussian product in 1D. ! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} END_DOC double precision , intent(in) :: a,b ! Exponents double precision , intent(in) :: xa,xb ! Centers double precision , intent(out) :: p ! New exponent double precision , intent(out) :: xp ! New center double precision , intent(out) :: k ! Constant double precision :: p_inv ASSERT (a>0.) ASSERT (b>0.) double precision :: xab, ab p = a+b p_inv = 1.d0/(a+b) ab = a*b xab = xa-xb ab = ab*p_inv k = ab*xab*xab if (k > 400.d0) then k=0.d0 return endif k = exp(-k) xp = (a*xa+b*xb)*p_inv end subroutine !- subroutine gaussian_product_x_v(a,xa,b,xb,k,p,xp,n_points) implicit none BEGIN_DOC ! Gaussian product in 1D with multiple xa ! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} END_DOC integer, intent(in) :: n_points double precision , intent(in) :: a,b ! Exponents double precision , intent(in) :: xa(n_points),xb ! Centers double precision , intent(out) :: p(n_points) ! New exponent double precision , intent(out) :: xp(n_points) ! New center double precision , intent(out) :: k(n_points) ! Constant double precision :: p_inv integer :: ipoint ASSERT (a>0.) ASSERT (b>0.) double precision :: xab, ab p = a+b p_inv = 1.d0/(a+b) ab = a*b*p_inv do ipoint = 1, n_points xab = xa(ipoint)-xb k(ipoint) = ab*xab*xab if (k(ipoint) > 40.d0) then k(ipoint)=0.d0 cycle endif k(ipoint) = exp(-k(ipoint)) xp(ipoint) = (a*xa(ipoint)+b*xb)*p_inv enddo end subroutine subroutine multiply_poly(b,nb,c,nc,d,nd) implicit none BEGIN_DOC ! Multiply two polynomials ! D(t) += B(t)*C(t) END_DOC integer, intent(in) :: nb, nc integer, intent(out) :: nd double precision, intent(in) :: b(0:nb), c(0:nc) double precision, intent(inout) :: d(0:nb+nc) integer :: ndtmp integer :: ib, ic, id, k if(ior(nc,nb) < 0) return !False if nc>=0 and nb>=0 select case (nb) case (0) call multiply_poly_b0(b,c,nc,d,nd) return case (1) call multiply_poly_b1(b,c,nc,d,nd) return case (2) call multiply_poly_b2(b,c,nc,d,nd) return end select select case (nc) case (0) call multiply_poly_c0(b,nb,c,d,nd) return case (1) call multiply_poly_c1(b,nb,c,d,nd) return case (2) call multiply_poly_c2(b,nb,c,d,nd) return end select do ib=0,nb do ic = 0,nc d(ib+ic) = d(ib+ic) + c(ic) * b(ib) enddo enddo do nd = nb+nc,0,-1 if (d(nd) /= 0.d0) exit enddo 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 ! 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 ! Add two polynomials ! D(t) += B(t)+C(t) END_DOC integer, intent(inout) :: nb, nc integer, intent(out) :: nd double precision, intent(in) :: b(0:nb), c(0:nc) double precision, intent(out) :: d(0:nb+nc) nd = nb+nc integer :: ib, ic, id do ib=0,max(nb,nc) d(ib) = d(ib) + c(ib) + b(ib) enddo do while ( (d(nd) == 0.d0).and.(nd>=0) ) nd -= 1 if (nd < 0) then exit endif enddo end subroutine add_poly_multiply(b,nb,cst,d,nd) implicit none BEGIN_DOC ! Add a polynomial multiplied by a constant ! D(t) += cst * B(t) END_DOC integer, intent(in) :: nb integer, intent(inout) :: nd double precision, intent(in) :: b(0:nb),cst double precision, intent(inout) :: d(0:max(nb,nd)) nd = max(nd,nb) if (nd /= -1) then integer :: ib, ic, id do ib=0,nb d(ib) = d(ib) + cst*b(ib) enddo do while ( d(nd) == 0.d0 ) nd -= 1 if (nd < 0) then exit endif enddo endif end subroutine recentered_poly2_v(P_new, lda, x_A, LD_xA, x_P, a, P_new2, ldb, x_B, x_Q, b, n_points) BEGIN_DOC ! Recenter two polynomials END_DOC implicit none integer, intent(in) :: a(3), b(3), n_points, lda, ldb, LD_xA double precision, intent(in) :: x_A(LD_xA,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(:,:) double precision :: fa, fb maxab(1:3) = max(a(1:3),b(1:3)) minab(1:3) = max(min(a(1:3),b(1:3)),(/0,0,0/)) allocate( pows_a(n_points,-2:maxval(maxab)+4), pows_b(n_points,-2:maxval(maxab)+4) ) do xyz=1,3 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(ipoint,xyz) - x_A(ipoint,xyz)) pows_b(ipoint,0) = 1.d0 pows_b(ipoint,1) = (x_Q(ipoint,xyz) - x_B(xyz)) enddo do i = 2,maxab(xyz) do ipoint=1,n_points pows_a(ipoint,i) = pows_a(ipoint,i-1)*pows_a(ipoint,1) pows_b(ipoint,i) = pows_b(ipoint,i-1)*pows_b(ipoint,1) enddo enddo do ipoint=1,n_points 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 (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 (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(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 (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(ipoint,i,xyz) = fb * pows_b(ipoint,b(xyz)-i) enddo enddo enddo 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 ! END_DOC implicit none integer, intent(in) :: a(3), n_points, lda, LD_xA double precision, intent(in) :: x_A(LD_xA,3), x_P(n_points,3) !double precision, intent(in) :: x_B(3), x_Q(n_points,3) double precision, intent(out) :: P_new(n_points,0:lda,3) !double precision, intent(out) :: P_new2(n_points,3) integer :: i, j, k, l, xyz, ipoint, maxab(3) double precision :: fa double precision, allocatable :: pows_a(:,:) !double precision, allocatable :: pows_b(:,:) double precision :: binom_func maxab(1:3) = max(a(1:3), (/0,0,0/)) allocate(pows_a(n_points,-2:maxval(maxab)+4)) !allocate(pows_b(n_points,-2:maxval(maxab)+4)) do xyz = 1, 3 if(a(xyz) < 0) cycle do ipoint = 1, n_points pows_a(ipoint,0) = 1.d0 pows_a(ipoint,1) = (x_P(ipoint,xyz) - x_A(ipoint,xyz)) !pows_b(ipoint,0) = 1.d0 !pows_b(ipoint,1) = (x_Q(ipoint,xyz) - x_B(xyz)) enddo do i = 2, maxab(xyz) do ipoint = 1, n_points pows_a(ipoint,i) = pows_a(ipoint,i-1) * pows_a(ipoint,1) !pows_b(ipoint,i) = pows_b(ipoint,i-1) * pows_b(ipoint,1) enddo enddo do ipoint = 1, n_points 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(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(ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i) enddo enddo enddo !xyz deallocate(pows_a) !deallocate(pows_b) end subroutine recentered_poly2_v0 subroutine recentered_poly2(P_new,x_A,x_P,a,P_new2,x_B,x_Q,b) implicit none BEGIN_DOC ! Recenter two polynomials END_DOC integer, intent(in) :: a,b double precision, intent(in) :: x_A,x_P,x_B,x_Q double precision, intent(out) :: P_new(0:a),P_new2(0:b) double precision :: pows_a(-2:a+b+4), pows_b(-2:a+b+4) double precision :: binom_func integer :: i,j,k,l, minab, maxab if ((a<0).or.(b<0) ) return maxab = max(a,b) minab = max(min(a,b),0) pows_a(0) = 1.d0 pows_a(1) = (x_P - x_A) pows_b(0) = 1.d0 pows_b(1) = (x_Q - x_B) do i = 2,maxab pows_a(i) = pows_a(i-1)*pows_a(1) pows_b(i) = pows_b(i-1)*pows_b(1) enddo P_new (0) = pows_a(a) P_new2(0) = pows_b(b) do i = 1,min(minab,20) P_new (i) = binom_transp(a-i,a) * pows_a(a-i) P_new2(i) = binom_transp(b-i,b) * pows_b(b-i) enddo do i = minab+1,min(a,20) P_new (i) = binom_transp(a-i,a) * pows_a(a-i) enddo do i = minab+1,min(b,20) P_new2(i) = binom_transp(b-i,b) * pows_b(b-i) enddo do i = 101,a P_new(i) = binom_func(a,a-i) * pows_a(a-i) enddo do i = 101,b P_new2(i) = binom_func(b,b-i) * pows_b(b-i) enddo end subroutine pol_modif_center(A_center, B_center, iorder, A_pol, B_pol) BEGIN_DOC ! ! Transform the pol centerd on A: ! [ \sum_i ax_i (x-x_A)^i ] [ \sum_j ay_j (y-y_A)^j ] [ \sum_k az_k (z-z_A)^k ] ! to a pol centered on B ! [ \sum_i bx_i (x-x_B)^i ] [ \sum_j by_j (y-y_B)^j ] [ \sum_k bz_k (z-z_B)^k ] ! END_DOC ! useful for max_dim include 'constants.include.F' implicit none integer, intent(in) :: iorder(3) double precision, intent(in) :: A_center(3), B_center(3) double precision, intent(in) :: A_pol(0:max_dim, 3) double precision, intent(out) :: B_pol(0:max_dim, 3) integer :: i, Lmax do i = 1, 3 Lmax = iorder(i) call pol_modif_center_x( A_center(i), B_center(i), Lmax, A_pol(0:Lmax, i), B_pol(0:Lmax, i) ) enddo return end subroutine pol_modif_center subroutine pol_modif_center_x(A_center, B_center, iorder, A_pol, B_pol) BEGIN_DOC ! ! Transform the pol centerd on A: ! [ \sum_i ax_i (x-x_A)^i ] ! to a pol centered on B ! [ \sum_i bx_i (x-x_B)^i ] ! ! bx_i = \sum_{j=i}^{iorder} ax_j (x_B - x_A)^(j-i) j! / [ i! (j-i)! ] ! = \sum_{j=i}^{iorder} ax_j (x_B - x_A)^(j-i) binom_func(j,i) ! END_DOC implicit none integer, intent(in) :: iorder double precision, intent(in) :: A_center, B_center double precision, intent(in) :: A_pol(0:iorder) double precision, intent(out) :: B_pol(0:iorder) integer :: i, j double precision :: fact_tmp, dx double precision :: binom_func dx = B_center - A_center do i = 0, iorder fact_tmp = 0.d0 do j = i, iorder fact_tmp += A_pol(j) * binom_func(j, i) * dx**dble(j-i) enddo B_pol(i) = fact_tmp enddo return end subroutine pol_modif_center_x double precision function F_integral(n,p) BEGIN_DOC ! function that calculates the following integral ! \int_{\-infty}^{+\infty} x^n \exp(-p x^2) dx END_DOC implicit none integer :: n double precision :: p integer :: i,j double precision :: accu,sqrt_p,fact_ratio,tmp,fact include 'constants.include.F' if(n < 0)then F_integral = 0.d0 endif if(iand(n,1).ne.0)then F_integral = 0.d0 return endif sqrt_p = 1.d0/dsqrt(p) if(n==0)then F_integral = sqpi * sqrt_p return endif F_integral = sqpi * 0.5d0**n * sqrt_p**(n+1) * fact(n)/fact(shiftr(n,1)) end double precision function rint(n,rho) implicit none BEGIN_DOC !.. math:: ! ! \int_0^1 dx \exp(-p x^2) x^n ! END_DOC include 'constants.include.F' double precision :: rho,u,rint1,v,val0,rint_large_n,u_inv integer :: n,k double precision :: two_rho_inv if(n.eq.0)then if(rho == 0.d0)then rint=1.d0 else u_inv=1.d0/dsqrt(rho) u=rho*u_inv rint=0.5d0*u_inv*sqpi*derf(u) endif return endif if(rho.lt.1.d0)then rint=rint1(n,rho) else if(n.le.20)then u_inv=1.d0/dsqrt(rho) if(rho.gt.80.d0)then v=0.d0 else v=dexp(-rho) endif u=rho*u_inv two_rho_inv = 0.5d0*u_inv*u_inv val0=0.5d0*u_inv*sqpi*derf(u) rint=(val0-v)*two_rho_inv do k=2,n rint=(rint*dfloat(k+k-1)-v)*two_rho_inv enddo else rint=rint_large_n(n,rho) endif endif end double precision function rint_sum(n_pt_out,rho,d1) implicit none BEGIN_DOC ! Needed for the calculation of two-electron integrals. END_DOC include 'constants.include.F' integer, intent(in) :: n_pt_out double precision, intent(in) :: rho,d1(0:n_pt_out) double precision :: u,rint1,v,val0,rint_large_n,u_inv integer :: n,k,i double precision :: two_rho_inv, rint_tmp, di if(rho < 1.d0)then if(rho == 0.d0)then rint_sum=d1(0) else u_inv=1.d0/dsqrt(rho) u=rho*u_inv rint_sum=0.5d0*u_inv*sqpi*derf(u) *d1(0) endif do i=2,n_pt_out,2 n = shiftr(i,1) rint_sum = rint_sum + d1(i)*rint1(n,rho) enddo else if(rho.gt.80.d0)then v=0.d0 else v=dexp(-rho) endif u_inv=1.d0/dsqrt(rho) u=rho*u_inv two_rho_inv = 0.5d0*u_inv*u_inv val0=0.5d0*u_inv*sqpi*derf(u) rint_sum=val0*d1(0) rint_tmp=(val0-v)*two_rho_inv di = 3.d0 do i=2,min(n_pt_out,40),2 rint_sum = rint_sum + d1(i)*rint_tmp rint_tmp = (rint_tmp*di-v)*two_rho_inv di = di+2.d0 enddo do i=42,n_pt_out,2 n = shiftr(i,1) rint_sum = rint_sum + d1(i)*rint_large_n(n,rho) enddo endif end double precision function hermite(n,x) implicit none BEGIN_DOC ! Hermite polynomial END_DOC integer :: n,k double precision :: h0,x,h1,h2 h0=1.d0 if(n.eq.0)then hermite=h0 return endif h1=x+x if(n.eq.1)then hermite=h1 return endif do k=1,n-1 h2=(x+x)*h1-dfloat(k+k)*h0 h0=h1 h1=h2 enddo hermite=h2 end double precision function rint_large_n(n,rho) implicit none BEGIN_DOC ! Version of rint for large values of n END_DOC integer :: n,k,l double precision :: rho,u,accu,eps,t1,t2,fact,alpha_k,rajout,hermite u=dsqrt(rho) accu=0.d0 k=0 eps=1.d0 do while (eps.gt.1.d-15) t1=1.d0 do l=0,k t1=t1*(n+n+l+1.d0) enddo t2=0.d0 do l=0,k t2=t2+(-1.d0)**l/(fact(l+1)*fact(k-l)) enddo alpha_k=t2*fact(k+1)*fact(k)*(-1.d0)**k alpha_k= alpha_k/t1 rajout=(-1.d0)**k*u**k*hermite(k,u)*alpha_k/fact(k) accu=accu+rajout eps=dabs(rajout)/accu k=k+1 enddo rint_large_n=dexp(-rho)*accu end double precision function rint1(n,rho) implicit none BEGIN_DOC ! Standard version of rint END_DOC integer, intent(in) :: n double precision, intent(in) :: rho double precision, parameter :: eps=1.d-15 double precision :: rho_tmp, diff integer :: k rint1=inv_int(n+n+1) rho_tmp = 1.d0 do k=1,20 rho_tmp = -rho_tmp*rho diff=rho_tmp*fact_inv(k)*inv_int(shiftl(k+n,1)+1) rint1=rint1+diff if (dabs(diff) > eps) then cycle endif return enddo write(*,*)'pb in rint1 k too large!' stop 1 end ! --- double precision function V_phi(n, m) BEGIN_DOC ! Computes the angular $\phi$ part of the nuclear attraction integral: ! ! $\int_{0}^{2 \pi} \cos(\phi)^n \sin(\phi)^m d\phi$. END_DOC implicit none integer, intent(in) :: n, m integer :: i double precision :: prod double precision :: Wallis prod = 1.d0 do i = 0, shiftr(n, 1)-1 prod = prod/ (1.d0 + dfloat(m+1)/dfloat(n-i-i-1)) enddo V_phi = 4.d0 * prod * Wallis(m) end function V_phi ! --- double precision function V_theta(n, m) BEGIN_DOC ! Computes the angular $\theta$ part of the nuclear attraction integral: ! ! $\int_{0}^{\pi} \cos(\theta)^n \sin(\theta)^m d\theta$ END_DOC implicit none include 'utils/constants.include.F' integer, intent(in) :: n, m integer :: i double precision :: prod double precision :: Wallis V_theta = 0.d0 prod = 1.d0 do i = 0, shiftr(n, 1)-1 prod = prod / (1.d0 + dfloat(m+1)/dfloat(n-i-i-1)) enddo V_theta = (prod + prod) * Wallis(m) end function V_theta ! --- double precision function Wallis(n) BEGIN_DOC ! Wallis integral: ! ! $\int_{0}^{\pi} \cos(\theta)^n d\theta$. END_DOC implicit none include 'utils/constants.include.F' integer, intent(in) :: n integer :: p double precision :: fact if(iand(n, 1) .eq. 0) then Wallis = fact(shiftr(n, 1)) Wallis = pi * fact(n) / (dble(ibset(0_8, n)) * (Wallis + Wallis) * Wallis) else p = shiftr(n, 1) Wallis = fact(p) Wallis = dble(ibset(0_8, p+p)) * Wallis * Wallis / fact(p+p+1) endif end function Wallis ! ---