mirror of
https://github.com/LCPQ/quantum_package
synced 2024-08-16 09:38:31 +02:00
Working Local
This commit is contained in:
parent
cca5ebc404
commit
64ea60c727
@ -42,9 +42,9 @@ then
|
|||||||
echo "ERROR"
|
echo "ERROR"
|
||||||
exit 1
|
exit 1
|
||||||
fi
|
fi
|
||||||
${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" $atoms
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" $atoms
|
||||||
|
|
||||||
|
cp He.dz_filipi.basis ${tmpfile}
|
||||||
|
echo ${tmpfile}
|
||||||
|
|
||||||
|
@ -550,7 +550,6 @@ double precision,allocatable :: array_R_loc(:,:,:)
|
|||||||
double precision,allocatable :: array_coefs(:,:,:,:,:,:)
|
double precision,allocatable :: array_coefs(:,:,:,:,:,:)
|
||||||
double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg
|
double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg
|
||||||
|
|
||||||
|
|
||||||
fourpi=4.d0*dacos(-1.d0)
|
fourpi=4.d0*dacos(-1.d0)
|
||||||
f=fourpi**1.5d0
|
f=fourpi**1.5d0
|
||||||
ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
|
ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
|
||||||
@ -567,6 +566,7 @@ double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg
|
|||||||
|
|
||||||
if(ac.eq.0.d0.and.bc.eq.0.d0)then
|
if(ac.eq.0.d0.and.bc.eq.0.d0)then
|
||||||
accu=0.d0
|
accu=0.d0
|
||||||
|
|
||||||
do k=1,klocmax
|
do k=1,klocmax
|
||||||
accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k))
|
accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k))
|
||||||
enddo
|
enddo
|
||||||
@ -1727,87 +1727,6 @@ end
|
|||||||
RETURN
|
RETURN
|
||||||
END
|
END
|
||||||
|
|
||||||
double precision FUNCTION GAMMLN(XX)
|
|
||||||
implicit double precision(a-h,o-z)
|
|
||||||
REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER
|
|
||||||
DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, &
|
|
||||||
-1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/
|
|
||||||
DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/
|
|
||||||
X=XX-ONE
|
|
||||||
TMP=X+FPF
|
|
||||||
TMP=(X+HALF)*DLOG(TMP)-TMP
|
|
||||||
SER=ONE
|
|
||||||
DO 11 J=1,6
|
|
||||||
X=X+ONE
|
|
||||||
SER=SER+COF(J)/X
|
|
||||||
11 CONTINUE
|
|
||||||
GAMMLN=TMP+DLOG(STP*SER)
|
|
||||||
RETURN
|
|
||||||
END
|
|
||||||
FUNCTION GAMMP(A,X)
|
|
||||||
implicit double precision(a-h,o-z)
|
|
||||||
IF(X.LT.0.d0.OR.A.LE.0.d0)PAUSE
|
|
||||||
IF(X.LT.A+1.d0)THEN
|
|
||||||
CALL GSER(GAMMP,A,X,GLN)
|
|
||||||
ELSE
|
|
||||||
CALL GCF(GAMMCF,A,X,GLN)
|
|
||||||
GAMMP=1.d0-GAMMCF
|
|
||||||
ENDIF
|
|
||||||
RETURN
|
|
||||||
END
|
|
||||||
SUBROUTINE GCF(GAMMCF,A,X,GLN)
|
|
||||||
implicit double precision(a-h,o-z)
|
|
||||||
PARAMETER (ITMAX=100,EPS=3.D-7)
|
|
||||||
GLN=GAMMLN(A)
|
|
||||||
GOLD=0.d0
|
|
||||||
A0=1.d0
|
|
||||||
A1=X
|
|
||||||
B0=0.d0
|
|
||||||
B1=1.d0
|
|
||||||
FAC=1.d0
|
|
||||||
DO 11 N=1,ITMAX
|
|
||||||
AN=DFLOAT(N)
|
|
||||||
ANA=AN-A
|
|
||||||
A0=(A1+A0*ANA)*FAC
|
|
||||||
B0=(B1+B0*ANA)*FAC
|
|
||||||
ANF=AN*FAC
|
|
||||||
A1=X*A0+ANF*A1
|
|
||||||
B1=X*B0+ANF*B1
|
|
||||||
IF(A1.NE.0.d0)THEN
|
|
||||||
FAC=1.d0/A1
|
|
||||||
G=B1*FAC
|
|
||||||
IF(DABS((G-GOLD)/G).LT.EPS)GO TO 1
|
|
||||||
GOLD=G
|
|
||||||
ENDIF
|
|
||||||
11 CONTINUE
|
|
||||||
PAUSE 'A TOO LARGE, ITMAX TOO SMALL'
|
|
||||||
1 GAMMCF=DEXP(-X+A*DLOG(X)-GLN)*G
|
|
||||||
RETURN
|
|
||||||
END
|
|
||||||
SUBROUTINE GSER(GAMSER,A,X,GLN)
|
|
||||||
implicit double precision(a-h,o-z)
|
|
||||||
PARAMETER (ITMAX=100,EPS=3.D-7)
|
|
||||||
GLN=GAMMLN(A)
|
|
||||||
IF(X.LE.0.d0)THEN
|
|
||||||
IF(X.LT.0.d0)PAUSE
|
|
||||||
GAMSER=0.d0
|
|
||||||
RETURN
|
|
||||||
ENDIF
|
|
||||||
AP=A
|
|
||||||
SUM=1.d0/A
|
|
||||||
DEL=SUM
|
|
||||||
DO 11 N=1,ITMAX
|
|
||||||
AP=AP+1.d0
|
|
||||||
DEL=DEL*X/AP
|
|
||||||
SUM=SUM+DEL
|
|
||||||
IF(DABS(DEL).LT.DABS(SUM)*EPS)GO TO 1
|
|
||||||
11 CONTINUE
|
|
||||||
PAUSE 'A TOO LARGE, ITMAX TOO SMALL'
|
|
||||||
1 GAMSER=SUM*DEXP(-X+A*DLOG(X)-GLN)
|
|
||||||
RETURN
|
|
||||||
END
|
|
||||||
|
|
||||||
|
|
||||||
double precision function coef_nk(n,k)
|
double precision function coef_nk(n,k)
|
||||||
implicit none
|
implicit none
|
||||||
integer n,k
|
integer n,k
|
||||||
|
@ -8,7 +8,7 @@
|
|||||||
double precision :: A_center(3),B_center(3),C_center(3)
|
double precision :: A_center(3),B_center(3),C_center(3)
|
||||||
integer :: power_A(3),power_B(3)
|
integer :: power_A(3),power_B(3)
|
||||||
integer :: i,j,k,l,n_pt_in,m
|
integer :: i,j,k,l,n_pt_in,m
|
||||||
double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult, Vloc
|
double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult, Vloc, Vpseudo
|
||||||
integer :: nucl_numC
|
integer :: nucl_numC
|
||||||
! Important for OpenMP
|
! Important for OpenMP
|
||||||
|
|
||||||
@ -38,13 +38,55 @@
|
|||||||
print*, "dz_k_ezfio", dz_k
|
print*, "dz_k_ezfio", dz_k
|
||||||
|
|
||||||
|
|
||||||
|
call ezfio_get_pseudo_v_k(v_k)
|
||||||
|
call ezfio_get_pseudo_n_k(n_k)
|
||||||
|
call ezfio_get_pseudo_dz_k(dz_k)
|
||||||
|
|
||||||
|
print*, "klocmax", klocmax
|
||||||
|
|
||||||
|
print*, "n_k_ezfio", n_k
|
||||||
|
print*, "v_k_ezfio",v_k
|
||||||
|
print*, "dz_k_ezfio", dz_k
|
||||||
|
|
||||||
|
|
||||||
|
!
|
||||||
|
! |\ | _ ._ | _ _ _. |
|
||||||
|
! | \| (_) | | | (_) (_ (_| |
|
||||||
|
!
|
||||||
|
|
||||||
|
!! Parameters of non local part of pseudo:
|
||||||
|
|
||||||
|
integer :: kmax,lmax
|
||||||
|
integer, allocatable :: n_kl(:,:)
|
||||||
|
double precision, allocatable :: v_kl(:,:), dz_kl(:,:)
|
||||||
|
|
||||||
|
call ezfio_get_pseudo_lmax(lmax)
|
||||||
|
call ezfio_get_pseudo_kmax(kmax)
|
||||||
|
lmax = lmax - 1
|
||||||
|
|
||||||
|
allocate(n_kl(kmax,0:lmax), v_kl(kmax,0:lmax), dz_kl(kmax,0:lmax))
|
||||||
|
|
||||||
|
call ezfio_get_pseudo_n_kl(n_kl)
|
||||||
|
call ezfio_get_pseudo_v_kl(v_kl)
|
||||||
|
call ezfio_get_pseudo_dz_kl(dz_kl)
|
||||||
|
|
||||||
|
|
||||||
|
print*, "lmax",lmax
|
||||||
|
print*, "kmax", kmax
|
||||||
|
|
||||||
|
print*,"n_kl_ezfio", n_kl
|
||||||
|
print*,"v_kl_ezfio", v_kl
|
||||||
|
print*,"dz_kl_ezfio", dz_kl
|
||||||
|
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP DEFAULT (NONE) &
|
||||||
!$OMP PRIVATE (i, j, k, l, m, alpha, beta, A_center, B_center, C_center, power_A, power_B, &
|
!$OMP PRIVATE (i, j, k, l, m, alpha, beta, A_center, B_center, C_center, power_A, power_B, &
|
||||||
!$OMP num_A, num_B, Z, c, n_pt_in,&
|
!$OMP num_A, num_B, Z, c, n_pt_in) &
|
||||||
!$OMP v_k, n_k, dz_k, klocmax) &
|
|
||||||
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, &
|
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, &
|
||||||
!$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge)
|
!$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge, &
|
||||||
|
!$OMP v_k, n_k, dz_k, klocmax, &
|
||||||
|
!$OMP lmax,kmax,v_kl,n_kl,dz_kl)
|
||||||
n_pt_in = n_pt_max_integrals
|
n_pt_in = n_pt_max_integrals
|
||||||
!$OMP DO SCHEDULE (guided)
|
!$OMP DO SCHEDULE (guided)
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
@ -74,15 +116,10 @@
|
|||||||
C_center(1) = nucl_coord(k,1)
|
C_center(1) = nucl_coord(k,1)
|
||||||
C_center(2) = nucl_coord(k,2)
|
C_center(2) = nucl_coord(k,2)
|
||||||
C_center(3) = nucl_coord(k,3)
|
C_center(3) = nucl_coord(k,3)
|
||||||
|
|
||||||
c = c + Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in)
|
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,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center)
|
||||||
print*, A_center
|
c = c - Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,C_center)
|
||||||
print*, B_center
|
|
||||||
print*, C_center
|
|
||||||
|
|
||||||
print*, Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center)
|
|
||||||
! c = c + Z*Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center)
|
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) - &
|
ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) - &
|
||||||
ao_coef_transp(l,j)*ao_coef_transp(m,i)*c
|
ao_coef_transp(l,j)*ao_coef_transp(m,i)*c
|
||||||
|
@ -1,289 +0,0 @@
|
|||||||
|
|
||||||
double precision function SABpartial(zA,zB,A,B,nA,nB,gamA,gamB)
|
|
||||||
implicit double precision(a-h,o-z)
|
|
||||||
dimension nA(3),nB(3)
|
|
||||||
dimension A(3),B(3)
|
|
||||||
gamtot=gamA+gamB
|
|
||||||
SABpartial=1.d0
|
|
||||||
|
|
||||||
l=3
|
|
||||||
u=gamA/gamtot*A(l)+gamB/gamtot*B(l)
|
|
||||||
arg=gamtot*u**2-gamA*A(l)**2-gamB*B(l)**2
|
|
||||||
alpha=dexp(arg)
|
|
||||||
&/gamtot**((1.d0+dfloat(nA(l))+dfloat(nB(l)))/2.d0)
|
|
||||||
wA=dsqrt(gamtot)*(u-A(l))
|
|
||||||
wB=dsqrt(gamtot)*(u-B(l))
|
|
||||||
boundA=dsqrt(gamtot)*(zA-u)
|
|
||||||
boundB=dsqrt(gamtot)*(zB-u)
|
|
||||||
|
|
||||||
accu=0.d0
|
|
||||||
do n=0,nA(l)
|
|
||||||
do m=0,nB(l)
|
|
||||||
integ=nA(l)+nB(l)-n-m
|
|
||||||
accu=accu
|
|
||||||
& +wA**n*wB**m*binom(nA(l),n)*binom(nB(l),m)
|
|
||||||
& *(rinteg(integ,boundB)-rinteg(integ,boundA))
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
SABpartial=SABpartial*accu*alpha
|
|
||||||
end
|
|
||||||
|
|
||||||
double precision function rintgauss(n)
|
|
||||||
implicit double precision(a-h,o-z)
|
|
||||||
rintgauss=dsqrt(dacos(-1.d0))
|
|
||||||
if(n.eq.0)return
|
|
||||||
if(n.eq.1)then
|
|
||||||
rintgauss=0.d0
|
|
||||||
return
|
|
||||||
endif
|
|
||||||
if(iand(n,1).eq.1)then
|
|
||||||
rintgauss=0.d0
|
|
||||||
return
|
|
||||||
endif
|
|
||||||
rintgauss=rintgauss/2.d0**(n/2)
|
|
||||||
rintgauss=rintgauss*ddfact2(n-1)
|
|
||||||
end
|
|
||||||
|
|
||||||
double precision function rinteg(n,u)
|
|
||||||
implicit double precision(a-h,o-z)
|
|
||||||
include 'constants.F'
|
|
||||||
! pi=dacos(-1.d0)
|
|
||||||
ichange=1
|
|
||||||
factor=1.d0
|
|
||||||
if(u.lt.0.d0)then
|
|
||||||
u=-u
|
|
||||||
factor=(-1.d0)**(n+1)
|
|
||||||
ichange=-1
|
|
||||||
endif
|
|
||||||
if(iand(n,1).eq.0)then
|
|
||||||
rinteg=0.d0
|
|
||||||
do l=0,n-2,2
|
|
||||||
prod=b_coef(l,u)
|
|
||||||
do k=l+2,n-2,2
|
|
||||||
prod=prod*a_coef(k)
|
|
||||||
enddo
|
|
||||||
rinteg=rinteg+prod
|
|
||||||
enddo
|
|
||||||
prod=dsqrt(pi)/2.d0*erf0(u)
|
|
||||||
do k=0,n-2,2
|
|
||||||
prod=prod*a_coef(k)
|
|
||||||
enddo
|
|
||||||
rinteg=rinteg+prod
|
|
||||||
endif
|
|
||||||
|
|
||||||
if(iand(n,1).eq.1)then
|
|
||||||
rinteg=0.d0
|
|
||||||
do l=1,n-2,2
|
|
||||||
prod=b_coef(l,u)
|
|
||||||
do k=l+2,n-2,2
|
|
||||||
prod=prod*a_coef(k)
|
|
||||||
enddo
|
|
||||||
rinteg=rinteg+prod
|
|
||||||
enddo
|
|
||||||
prod=0.5d0*(1.d0-dexp(-u**2))
|
|
||||||
do k=1,n-2,2
|
|
||||||
prod=prod*a_coef(k)
|
|
||||||
enddo
|
|
||||||
rinteg=rinteg+prod
|
|
||||||
endif
|
|
||||||
|
|
||||||
rinteg=rinteg*factor
|
|
||||||
|
|
||||||
if(ichange.eq.-1)u=-u
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
!<function type="double precision function" name="erf0">
|
|
||||||
! <arg name="x"
|
|
||||||
! doc ="" />
|
|
||||||
!
|
|
||||||
! <doc>
|
|
||||||
!
|
|
||||||
! </doc>
|
|
||||||
!
|
|
||||||
! <fortran>
|
|
||||||
double precision function erf0(x)
|
|
||||||
implicit double precision (a-h,o-z)
|
|
||||||
if(x.lt.0.d0)then
|
|
||||||
erf0=-gammp(0.5d0,x**2)
|
|
||||||
else
|
|
||||||
erf0=gammp(0.5d0,x**2)
|
|
||||||
endif
|
|
||||||
end
|
|
||||||
|
|
||||||
|
|
||||||
! </fortran>
|
|
||||||
!</function>
|
|
||||||
!<function type="double precision function" name="gammp">
|
|
||||||
! <arg name="a"
|
|
||||||
! doc ="" />
|
|
||||||
! <arg name="x"
|
|
||||||
! doc ="" />
|
|
||||||
!
|
|
||||||
! <doc>
|
|
||||||
!
|
|
||||||
! </doc>
|
|
||||||
!
|
|
||||||
! <calls>
|
|
||||||
! gcf
|
|
||||||
! gser
|
|
||||||
! </calls>
|
|
||||||
!
|
|
||||||
! <fortran>
|
|
||||||
double precision function gammp(a,x)
|
|
||||||
implicit double precision (a-h,o-z)
|
|
||||||
if(x.lt.0..or.a.le.0.)stop 'error in gammp'
|
|
||||||
if(x.lt.a+1.)then
|
|
||||||
call gser(gammp,a,x,gln)
|
|
||||||
else
|
|
||||||
call gcf(gammcf,a,x,gln)
|
|
||||||
gammp=1.-gammcf
|
|
||||||
endif
|
|
||||||
return
|
|
||||||
end
|
|
||||||
! </fortran>
|
|
||||||
!</function>
|
|
||||||
|
|
||||||
|
|
||||||
!<function type="subroutine" name="gser">
|
|
||||||
! <arg name="gamser"
|
|
||||||
! doc ="" />
|
|
||||||
! <arg name="a"
|
|
||||||
! doc ="" />
|
|
||||||
! <arg name="x"
|
|
||||||
! doc ="" />
|
|
||||||
! <arg name="gln"
|
|
||||||
! doc ="" />
|
|
||||||
!
|
|
||||||
! <doc>
|
|
||||||
!
|
|
||||||
! </doc>
|
|
||||||
!
|
|
||||||
! <calledBy>
|
|
||||||
! gammp
|
|
||||||
! </calledBy>
|
|
||||||
!
|
|
||||||
! <fortran>
|
|
||||||
subroutine gser(gamser,a,x,gln)
|
|
||||||
implicit double precision (a-h,o-z)
|
|
||||||
parameter (itmax=100,eps=3.e-7)
|
|
||||||
gln=gammln(a)
|
|
||||||
if(x.le.0.)then
|
|
||||||
if(x.lt.0.) stop 'error in gser'
|
|
||||||
gamser=0.
|
|
||||||
return
|
|
||||||
endif
|
|
||||||
ap=a
|
|
||||||
sum=1./a
|
|
||||||
del=sum
|
|
||||||
do 11 n=1,itmax
|
|
||||||
ap=ap+1.
|
|
||||||
del=del*x/ap
|
|
||||||
sum=sum+del
|
|
||||||
if(abs(del).lt.abs(sum)*eps)go to 1
|
|
||||||
11 continue
|
|
||||||
stop 'a too large, itmax too small'
|
|
||||||
1 gamser=sum*exp(-x+a*log(x)-gln)
|
|
||||||
return
|
|
||||||
end
|
|
||||||
! </fortran>
|
|
||||||
|
|
||||||
!</function>
|
|
||||||
!<function type="subroutine" name="gcf">
|
|
||||||
! <arg name="gammcf"
|
|
||||||
! doc ="" />
|
|
||||||
! <arg name="a"
|
|
||||||
! doc ="" />
|
|
||||||
! <arg name="x"
|
|
||||||
! doc ="" />
|
|
||||||
! <arg name="gln"
|
|
||||||
! doc ="" />
|
|
||||||
!
|
|
||||||
! <doc>
|
|
||||||
!
|
|
||||||
! </doc>
|
|
||||||
!
|
|
||||||
! <calledBy>
|
|
||||||
! gammp
|
|
||||||
! </calledBy>
|
|
||||||
!
|
|
||||||
! <fortran>
|
|
||||||
subroutine gcf(gammcf,a,x,gln)
|
|
||||||
implicit double precision (a-h,o-z)
|
|
||||||
parameter (itmax=100,eps=3.e-7)
|
|
||||||
gln=gammln(a)
|
|
||||||
gold=0.
|
|
||||||
a0=1.
|
|
||||||
a1=x
|
|
||||||
b0=0.
|
|
||||||
b1=1.
|
|
||||||
fac=1.
|
|
||||||
do 11 n=1,itmax
|
|
||||||
an=float(n)
|
|
||||||
ana=an-a
|
|
||||||
a0=(a1+a0*ana)*fac
|
|
||||||
b0=(b1+b0*ana)*fac
|
|
||||||
anf=an*fac
|
|
||||||
a1=x*a0+anf*a1
|
|
||||||
b1=x*b0+anf*b1
|
|
||||||
if(a1.ne.0.)then
|
|
||||||
fac=1./a1
|
|
||||||
g=b1*fac
|
|
||||||
if(abs((g-gold)/g).lt.eps)go to 1
|
|
||||||
gold=g
|
|
||||||
endif
|
|
||||||
11 continue
|
|
||||||
stop 'a too large, itmax too small'
|
|
||||||
1 gammcf=exp(-x+a*log(x)-gln)*g
|
|
||||||
return
|
|
||||||
end
|
|
||||||
|
|
||||||
! </fortran>
|
|
||||||
!</function>
|
|
||||||
double precision function ddfact2(n)
|
|
||||||
implicit double precision(a-h,o-z)
|
|
||||||
if(iand(n,1).eq.0)stop 'error in ddfact2'
|
|
||||||
ddfact2=1.d0
|
|
||||||
do i=1,n,2
|
|
||||||
ddfact2=ddfact2*dfloat(i)
|
|
||||||
enddo
|
|
||||||
end
|
|
||||||
|
|
||||||
double precision function a_coef(n)
|
|
||||||
implicit double precision(a-h,o-z)
|
|
||||||
a_coef=dfloat(n+1)/2.d0
|
|
||||||
end
|
|
||||||
|
|
||||||
double precision function b_coef(n,u)
|
|
||||||
implicit double precision(a-h,o-z)
|
|
||||||
b_coef=-0.5d0*u**(n+1)*dexp(-u**2)
|
|
||||||
end
|
|
||||||
|
|
||||||
!<function type="double precision function" name="gammln">
|
|
||||||
! <arg name="xx"
|
|
||||||
! doc ="" />
|
|
||||||
!
|
|
||||||
! <doc>
|
|
||||||
!
|
|
||||||
! </doc>
|
|
||||||
!
|
|
||||||
! <fortran>
|
|
||||||
double precision function gammln(xx)
|
|
||||||
implicit double precision (a-h,o-z)
|
|
||||||
real*8 cof(6),stp,half,one,fpf,x,tmp,ser
|
|
||||||
data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
|
|
||||||
* -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
|
|
||||||
data half,one,fpf/0.5d0,1.0d0,5.5d0/
|
|
||||||
x=xx-one
|
|
||||||
tmp=x+fpf
|
|
||||||
tmp=(x+half)*log(tmp)-tmp
|
|
||||||
ser=one
|
|
||||||
do 11 j=1,6
|
|
||||||
x=x+one
|
|
||||||
ser=ser+cof(j)/x
|
|
||||||
11 continue
|
|
||||||
gammln=tmp+log(stp*ser)
|
|
||||||
return
|
|
||||||
end
|
|
||||||
! </fortran>
|
|
||||||
!</function>
|
|
Loading…
Reference in New Issue
Block a user