10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-07-11 22:03:47 +02:00

>50% acceleration of int.f90

This commit is contained in:
Thomas Applencourt 2015-05-05 14:22:38 +02:00
parent db8f905b98
commit c1c684b0e7
2 changed files with 129 additions and 96 deletions

View File

@ -1898,114 +1898,147 @@ end
!! !!
!! M_n(x) modified spherical bessel function !! M_n(x) modified spherical bessel function
!! !!
double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
implicit none
integer n,k,m,q,l,kcp
double precision gam,dblefact,fact,pi,a,b
double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
logical done
u=(a+b)/(2.d0*dsqrt(gam)) double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
freal=dexp(-arg)
if(a.eq.0.d0.and.b.eq.0.d0)then implicit none
if(n.ne.0.or.m.ne.0)then integer n,k,m,q,l,kcp
int_prod_bessel=0.d0 double precision gam,dblefact,fact,pi,a,b
return double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
endif double precision dump
int_prod_bessel=crochet(l,gam)*freal
return
endif
if(u.gt.6.d0)then double precision, allocatable :: temp_array(:)
int_prod_bessel=int_prod_bessel_large(l,gam,n,m,a,b,arg)
logical done
u=(a+b)/(2.d0*dsqrt(gam))
freal=dexp(-arg)
if(a.eq.0.d0.and.b.eq.0.d0)then
if(n.ne.0.or.m.ne.0)then
int_prod_bessel=0.d0
return return
endif
int_prod_bessel=crochet(l,gam)*freal
return
endif
if(u.gt.6.d0)then
int_prod_bessel=int_prod_bessel_large(l,gam,n,m,a,b,arg)
return
endif
if(a.ne.0.d0.and.b.ne.0.d0)then
q=0
intold=-1.d0
int=0.d0
done=.false.
kcp=0
allocate( temp_array(0:101))
do while (.not.done)
kcp=kcp+1
sum=0.d0
dump = a**n
do k=q,q+1
temp_array(k) = coef_nk(n,k)*dump*(a**(2*k))
enddo
do k=0,q
sum=sum+temp_array(k)*coef_nk(m,q-k)*b**(m-2*k+2*q)
enddo
int=int+sum*crochet(2*q+n+m+l,gam)
if(dabs(int-intold).lt.1d-15)then
done=.true.
else
q=q+1
intold=int
endif endif
if(a.ne.0.d0.and.b.ne.0.d0)then enddo
deallocate(temp_array)
q=0 int_prod_bessel=int*freal
intold=-1.d0
int=0.d0 if(kcp.gt.100)then
done=.false. print*,'**WARNING** bad convergence in int_prod_bessel'
kcp=0 endif
do while (.not.done)
kcp=kcp+1 return
sum=0.d0 endif
do k=0,q
sum=sum+coef_nk(n,k)*coef_nk(m,q-k)*a**(n+2*k)*b**(m-2*k+2*q) if(a.eq.0.d0.and.b.ne.0.d0)then
enddo
int=int+sum*crochet(2*q+n+m+l,gam) if(n.ne.0)then
if(dabs(int-intold).lt.1d-15)then int_prod_bessel=0.d0
return
endif
q=0
intold=-1.d0
int=0.d0
done=.false.
kcp=0
do while (.not.done)
kcp=kcp+1
int=int+coef_nk(m,q)*b**(m+2*q)*crochet(2*q+m+l,gam)
if(dabs(int-intold).lt.1d-15)then
done=.true.
else
q=q+1
intold=int
endif
enddo
int_prod_bessel=int*freal
if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
return
endif
if(a.ne.0.d0.and.b.eq.0.d0)then
if(m.ne.0)then
int_prod_bessel=0.d0
return
endif
q=0
intold=-1.d0
int=0.d0
done=.false.
kcp=0
do while (.not.done)
kcp=kcp+1
int=int+coef_nk(n,q)*a**(n+2*q)*crochet(2*q+n+l,gam)
if(dabs(int-intold).lt.1d-15)then
done=.true. done=.true.
else else
q=q+1 q=q+1
intold=int intold=int
endif
enddo
int_prod_bessel=int*freal
if(kcp.gt.100)then
print*,'**WARNING** bad convergence in int_prod_bessel'
! else
! print*,'kcp=',kcp
endif
return
endif endif
enddo
int_prod_bessel=int*freal
if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
return
if(a.eq.0.d0.and.b.ne.0.d0)then endif
if(n.ne.0)then
int_prod_bessel=0.d0
return
endif
q=0 stop 'pb in int_prod_bessel!!'
intold=-1.d0 end
int=0.d0
done=.false.
kcp=0
do while (.not.done)
kcp=kcp+1
int=int+coef_nk(m,q)*b**(m+2*q)*crochet(2*q+m+l,gam)
if(dabs(int-intold).lt.1d-15)then
done=.true.
else
q=q+1
intold=int
endif
enddo
int_prod_bessel=int*freal
if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
return
endif
if(a.ne.0.d0.and.b.eq.0.d0)then
if(m.ne.0)then
int_prod_bessel=0.d0
return
endif
q=0
intold=-1.d0
int=0.d0
done=.false.
kcp=0
do while (.not.done)
kcp=kcp+1
int=int+coef_nk(n,q)*a**(n+2*q)*crochet(2*q+n+l,gam)
if(dabs(int-intold).lt.1d-15)then
done=.true.
else
q=q+1
intold=int
endif
enddo
int_prod_bessel=int*freal
if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
return
endif
stop 'pb in int_prod_bessel!!'
end
double precision function int_prod_bessel_large(l,gam,n,m,a,b,arg) double precision function int_prod_bessel_large(l,gam,n,m,a,b,arg)
implicit none implicit none

View File

@ -169,7 +169,7 @@ def run_hf(geo, basis, mult=1, pseudo=False, remove_after_sucess=True):
ref_energy["sto-3g"]["methane"] = Energy(-39.7267433402, None) ref_energy["sto-3g"]["methane"] = Energy(-39.7267433402, None)
ref_energy["vdz"]["SO2"] = Energy(None, -41.48912297776174) ref_energy["vdz"]["SO2"] = Energy(None, -41.48912297776174)
ref_energy["vdz"]["HBO"] = Energy(None, -19.1198251160) ref_energy["vdz"]["HBO"] = Energy(None, -19.11982540413317)
# ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
# G l o b a l _ v a r i a b l e # # G l o b a l _ v a r i a b l e #