10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-26 07:02:14 +02:00

Merge branch 'master' of github.com:scemama/quantum_package

This commit is contained in:
Anthony Scemama 2016-09-30 15:35:10 +02:00
commit b21cfe99c8
3 changed files with 29 additions and 19 deletions

View File

@ -21,9 +21,9 @@ use bitmasks
do k=1,N_int
good = good .and. ( &
iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == &
iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( &
iand(not(cas_bitmask(k,1,l)), hf_bitmask(k,1)) ) .and. ( &
iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == &
iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) )
iand(not(cas_bitmask(k,2,l)), hf_bitmask(k,2)) )
enddo
if (good) then
exit

View File

@ -136,7 +136,7 @@ end
if(l.eq.2.and.m.eq.2)ylm_real=dsqrt(15.d0/16.d0/pi)*(xchap*xchap-ychap*ychap)
if(l.eq.2.and.m.eq.1)ylm_real=dsqrt(15.d0/fourpi)*xchap*zchap
if(l.eq.2.and.m.eq.0)ylm_real=dsqrt(5.d0/16.d0/pi)*(-xchap*xchap-ychap*ychap+2.d0*zchap*zchap)
if(l.eq.2.and.m.eq.0)ylm_real=dsqrt(5.d0/16.d0/pi)*(2.d0*zchap*zchap-xchap*xchap-ychap*ychap)
if(l.eq.2.and.m.eq.-1)ylm_real=dsqrt(15.d0/fourpi)*ychap*zchap
if(l.eq.2.and.m.eq.-2)ylm_real=dsqrt(15.d0/fourpi)*xchap*ychap
@ -279,8 +279,7 @@ 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)*int_prod_bessel(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
@ -1783,10 +1782,8 @@ double precision function coef_nk(n,k)
double precision gam,dble_fact,fact
gam=dble_fact(n+n+k+k+1)
! coef_nk=1.d0/(dble(ISHFT(1,k))*fact(k)*gam)
coef_nk=1.d0/(2.d0**k*fact(k)*gam)
! coef_nk=1.d0/(2.d0**k*fact(k)*gam)
coef_nk=1.d0/(dble(ibset(0,k))*fact(k)*gam)
return
@ -1815,7 +1812,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
logical done
u=(a+b)/(2.d0*dsqrt(gam))
u=(a+b)*0.5d0/dsqrt(gam)
freal=dexp(-arg)
if(a.eq.0.d0.and.b.eq.0.d0)then

View File

@ -84,10 +84,8 @@ double precision function fact(n)
memo(i) = memo(i-1)*dble(i)
enddo
memomax = min(n,100)
fact = memo(memomax)
do i=101,n
fact = fact*dble(i)
enddo
double precision :: logfact
fact = dexp(logfact(n))
end function
double precision function logfact(n)
@ -164,12 +162,27 @@ double precision function dble_fact_even(n) result(fact2)
ASSERT (iand(n,1) /= 1)
prod=1.d0
do k=2,n,2
prod=prod*dfloat(k)
if (n <= memomax) then
if (n < 2) then
fact2 = 1.d0
else
fact2 = memo(n)
endif
return
endif
integer :: i
memo(1)=1.d0
do i=memomax+2,min(n,100),2
memo(i) = memo(i-2)* dble(i)
enddo
fact2=prod
return
memomax = min(n,100)
fact2 = memo(memomax)
if (n > 100) then
double precision :: dble_logfact
fact2 = dexp(dble_logfact(n))
endif
end function