diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index c7d2ac84..750dbaa7 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -140,7 +140,7 @@ end ! __ __ _ __ ___ ___ _ _ __| | ___ ! \ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ ! \ V / | |_) \__ \ __/ |_| | (_| | (_) | -! \_/ | .__/|___/\___|\__,_|\__,_|\___/ +! \_/ | .__/|___/\___|\__,_|\____|\___/ ! | | ! |_| @@ -200,7 +200,7 @@ double precision, intent(in) :: v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_ double precision :: fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big -double precision :: areal,freal,breal,t1,t2,int_prod_bessel +double precision :: areal,freal,breal,t1,t2,int_prod_bessel, int_prod_bessel_num_soph_p double precision :: arg integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot @@ -237,7 +237,7 @@ fourpi=4.d0*dacos(-1.d0) ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) arg=g_a*ac**2+g_b*bc**2 -if(arg.gt.-dlog(10.d-20))then +if(arg.gt.-dlog(1.d-20))then Vpseudo=0.d0 return endif @@ -270,7 +270,9 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then do m=-l,l prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) - accu=accu+prod*prodp*v_kl(k,l)*freal*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal) + + accu=accu+prod*prodp*v_kl(k,l)*int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg) + enddo enddo enddo @@ -302,8 +304,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then do lambdap=0,lmax+ntotB do k=1,kmax do l=0,lmax - array_R(ktot,k,l,lambda,lambdap)= freal & - *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal) + array_R(ktot,k,l,lambda,lambdap)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg) enddo enddo enddo @@ -425,9 +426,8 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then do k=1,kmax do l=0,lmax - array_R(ktot,k,l,0,lambdap)= freal & - *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal) - enddo + array_R(ktot,k,l,0,lambdap)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg) + enddo enddo enddo enddo @@ -512,9 +512,7 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then do k=1,kmax do l=0,lmax - array_R(ktot,k,l,lambda,0)= freal & - *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal) - + array_R(ktot,k,l,lambda,0)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg) enddo enddo enddo @@ -1974,6 +1972,61 @@ end stop 'pb in int_prod_bessel!!' end + +double precision function int_prod_bessel_num_soph_p(l,gam,n,m,a,b,arg) + implicit none + integer :: n,m,l + double precision :: gam,a,b,arg,arg_new + double precision :: bessel_mod,factor + logical :: not_done + double precision :: bigA,xold,x,dx,accu,intnew,intold,intold2,u,v,freal + integer :: iter, i, nI, n0 + double precision :: eps + + u=(a+b)/(2.d0*dsqrt(gam)) + arg_new=u**2-arg + freal=dexp(arg_new) + v=u/dsqrt(gam) + + bigA=v+dsqrt(-dlog(1.d-15)/gam) + n0=5 + accu=0.d0 + dx=bigA/(float(n0)-1.d0) + iter=0 + do i=1,n0 + x=(float(i)-1.d0)*dx + accu=accu+x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x) + enddo + +accu=accu*freal +intold=accu*dx + +eps=1.d-08 +nI=n0-1 +dx=dx/2.d0 +not_done=.true. + +do while(not_done) + iter=iter+1 + accu=0.d0 + do i=1,nI + x=dx+(float(i)-1.d0)*2.d0*dx + accu=accu+dx*x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x) + enddo + accu=accu*freal + intnew=intold/2.d0+accu + if(iter.gt.1.and.dabs(intnew-intold).lt.eps.and.dabs(intnew-intold2).lt.eps)then + not_done=.false. + else + intold2=intold + intold=intnew + dx=dx/2.d0 + nI=2*nI + endif +enddo +int_prod_bessel_num_soph_p=intold +end + !! Calculation of !! !! I= \int dx x**l *exp(-gam*x**2) M_n(ax) diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index ef2e8c8f..badf4afd 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -156,7 +156,7 @@ c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) - c = c + Vloc( klocmax ,v_k(k,:) ,n_k(k,:) ,dz_k(k,:), A_center,power_A,alpha,B_center,power_B,beta,C_center) +! c = c + Vloc( klocmax ,v_k(k,:) ,n_k(k,:) ,dz_k(k,:), A_center,power_A,alpha,B_center,power_B,beta,C_center) n_kl_dump = n_kl(k,1:kmax,0:lmax) @@ -167,7 +167,7 @@ ! print*, "kmax",kmax ! print*, "v_kl",v_kl_dump ! print*, "n_kl",n_kl_dump -! print*, n_kl_dump(1,0) +! print*, n_kl_ump(1,0) ! print*, n_kl_dump(1,1) ! print*, "dz_kl",dz_kl_dump ! print*, dz_kl_dump(1,0) @@ -180,7 +180,9 @@ ! print*, "beta", beta ! print*, "C_center",C_center - ! c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) + c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) + dump = Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) + print*, dump ! c = c - Vps(A_center,power_A,alpha,B_center,power_B,beta,C_center,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) ! print*, "#################"