LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-04 05:03:54 +01:00
Code Issues Releases Wiki Activity

HF_bitmask in cas definition

Browse Source
This commit is contained in:
Anthony Scemama 2016-09-30 15:34:59 +02:00
parent 7219dda33b
commit 41e69c56e5
3 changed files with 29 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/Determinants/psi_cas.irp.f
View File

@ -21,9 +21,9 @@ use bitmasks
do k=1,N_int do k=1,N_int
good = good .and. ( & 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,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,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 enddo
if (good) then if (good) then
exit exit

13
src/Integrals_Monoelec/pseudopot.f90
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.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.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.-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 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 do m=-l,l
prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) 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)) 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
enddo enddo
@ -1783,10 +1782,8 @@ double precision function coef_nk(n,k)
double precision gam,dble_fact,fact double precision gam,dble_fact,fact
gam=dble_fact(n+n+k+k+1) gam=dble_fact(n+n+k+k+1)
! coef_nk=1.d0/(2.d0**k*fact(k)*gam)
! coef_nk=1.d0/(dble(ISHFT(1,k))*fact(k)*gam) coef_nk=1.d0/(dble(ibset(0,k))*fact(k)*gam)
coef_nk=1.d0/(2.d0**k*fact(k)*gam)
return return
@ -1815,7 +1812,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
logical done logical done
u=(a+b)/(2.d0*dsqrt(gam)) u=(a+b)*0.5d0/dsqrt(gam)
freal=dexp(-arg) freal=dexp(-arg)
if(a.eq.0.d0.and.b.eq.0.d0)then if(a.eq.0.d0.and.b.eq.0.d0)then

31
src/Utils/util.irp.f
View File

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