mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-23 04:43:50 +01:00
Integration with the new code
This commit is contained in:
Thomas Applencourt
parent
eeb60de46d
commit
db8f905b98
@ -195,7 +195,7 @@ double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,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, int_prod_bessel_num_soph_p
|
||||
double precision :: areal,freal,breal,t1,t2,int_prod_bessel
|
||||
double precision :: arg
|
||||
|
||||
integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot
|
||||
@ -276,7 +276,7 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then
|
||||
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)*int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg)
|
||||
accu=accu+prod*prodp*v_kl(k,l)*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
@ -309,7 +309,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)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg)
|
||||
array_R(ktot,k,l,lambda,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -431,7 +431,7 @@ 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)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
|
||||
array_R(ktot,k,l,0,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -517,7 +517,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)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg)
|
||||
array_R(ktot,k,l,lambda,0)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -1898,14 +1898,32 @@ end
|
||||
!!
|
||||
!! M_n(x) modified spherical bessel function
|
||||
!!
|
||||
double precision function int_prod_bessel(l,gam,n,m,a,b)
|
||||
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
|
||||
double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,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
|
||||
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
|
||||
@ -1925,17 +1943,12 @@ end
|
||||
intold=int
|
||||
endif
|
||||
enddo
|
||||
int_prod_bessel=int
|
||||
if(kcp.gt.100) then
|
||||
print*,"l",l
|
||||
print*, "gam", gam
|
||||
print*, "n", n
|
||||
print*, "m", m
|
||||
print*, "a", a
|
||||
print*, "b", b
|
||||
print*, "kcp", kcp
|
||||
print*,'**WARNING** bad convergence in int_prod_bessel'
|
||||
endif
|
||||
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
|
||||
|
||||
@ -1944,6 +1957,7 @@ end
|
||||
int_prod_bessel=0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
q=0
|
||||
intold=-1.d0
|
||||
int=0.d0
|
||||
@ -1959,7 +1973,7 @@ end
|
||||
intold=int
|
||||
endif
|
||||
enddo
|
||||
int_prod_bessel=int
|
||||
int_prod_bessel=int*freal
|
||||
if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
|
||||
return
|
||||
endif
|
||||
@ -1969,6 +1983,7 @@ end
|
||||
int_prod_bessel=0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
q=0
|
||||
intold=-1.d0
|
||||
int=0.d0
|
||||
@ -1984,76 +1999,74 @@ end
|
||||
intold=int
|
||||
endif
|
||||
enddo
|
||||
int_prod_bessel=int
|
||||
int_prod_bessel=int*freal
|
||||
if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
|
||||
return
|
||||
endif
|
||||
|
||||
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
|
||||
endif
|
||||
int_prod_bessel=crochet(l,gam)
|
||||
return
|
||||
endif
|
||||
|
||||
stop 'pb in int_prod_bessel!!'
|
||||
end
|
||||
|
||||
|
||||
double precision function int_prod_bessel_num_soph_p(l,gam,n,m,a,b,arg)
|
||||
double precision function int_prod_bessel_large(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
|
||||
integer n,m,i,npts,l
|
||||
double precision gam,a,b
|
||||
double precision sum,x,bessel_mod,u,factor,arg
|
||||
double precision xq(100),wq(100)
|
||||
|
||||
u=(a+b)/(2.d0*dsqrt(gam))
|
||||
arg_new=u**2-arg
|
||||
freal=dexp(arg_new)
|
||||
v=u/dsqrt(gam)
|
||||
factor=dexp(u**2-arg)/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
|
||||
xq(1)= 5.38748089001123
|
||||
xq(2)= 4.60368244955074
|
||||
xq(3)= 3.94476404011563
|
||||
xq(4)= 3.34785456738322
|
||||
xq(5)= 2.78880605842813
|
||||
xq(6)= 2.25497400208928
|
||||
xq(7)= 1.73853771211659
|
||||
xq(8)= 1.23407621539532
|
||||
xq(9)= 0.737473728545394
|
||||
xq(10)= 0.245340708300901
|
||||
xq(11)=-0.245340708300901
|
||||
xq(12)=-0.737473728545394
|
||||
xq(13)=-1.23407621539532
|
||||
xq(14)=-1.73853771211659
|
||||
xq(15)=-2.25497400208928
|
||||
xq(16)=-2.78880605842813
|
||||
xq(17)=-3.34785456738322
|
||||
xq(18)=-3.94476404011563
|
||||
xq(19)=-4.60368244955074
|
||||
xq(20)=-5.38748089001123
|
||||
wq(1)= 2.229393645534151E-013
|
||||
wq(2)= 4.399340992273176E-010
|
||||
wq(3)= 1.086069370769280E-007
|
||||
wq(4)= 7.802556478532063E-006
|
||||
wq(5)= 2.283386360163528E-004
|
||||
wq(6)= 3.243773342237853E-003
|
||||
wq(7)= 2.481052088746362E-002
|
||||
wq(8)= 0.109017206020022
|
||||
wq(9)= 0.286675505362834
|
||||
wq(10)= 0.462243669600610
|
||||
wq(11)= 0.462243669600610
|
||||
wq(12)= 0.286675505362834
|
||||
wq(13)= 0.109017206020022
|
||||
wq(14)= 2.481052088746362E-002
|
||||
wq(15)= 3.243773342237853E-003
|
||||
wq(16)= 2.283386360163528E-004
|
||||
wq(17)= 7.802556478532063E-006
|
||||
wq(18)= 1.086069370769280E-007
|
||||
wq(19)= 4.399340992273176E-010
|
||||
wq(20)= 2.229393645534151E-013
|
||||
|
||||
accu=accu*freal
|
||||
intold=accu*dx
|
||||
npts=20
|
||||
! call gauher(xq,wq,npts)
|
||||
|
||||
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
|
||||
sum=0.d0
|
||||
do i=1,npts
|
||||
x=(xq(i)+u)/dsqrt(gam)
|
||||
sum=sum+wq(i)*x**l*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x)
|
||||
enddo
|
||||
int_prod_bessel_large=sum*factor
|
||||
end
|
||||
|
||||
!! Calculation of
|
||||
|
||||
@ -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["vdz"]["SO2"] = Energy(None, -41.48912297776174)
|
||||
ref_energy["vdz"]["HBO"] = Energy(None, -19.11982537379352)
|
||||
ref_energy["vdz"]["HBO"] = Energy(None, -19.1198251160)
|
||||
|
||||
# ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
|
||||
# G l o b a l _ v a r i a b l e #
|
||||
|
||||
Loading…
Reference in New Issue
Block a user