LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-10 21:53:29 +01:00
Code Releases Activity
409aa82373
qp2/src/ao_one_e_ints/pseudopot.f90

2128 lines
52 KiB
Fortran
Raw Normal View History

Initial commit
2019-01-25 11:39:31 +01:00
!! INFO : You can display equations using : http://www.codecogs.com/latex/eqneditor.php
!!
!! {\tt Vps}(C) = \langle \Phi_A|{\tt Vloc}(C)+{\tt Vpp}(C)| \Phi_B \rangle
!!
!! with :
!!
!! {\tt Vloc}(C)=\sum_{k=1}^{\tt klocmax} v_k r_C^{n_k} \exp(-dz_k r_C^2) \\
!!
!! {\tt Vpp}(C)=\sum_{l=0}^{\tt lmax}\left( \sum_{k=1}^{\tt kmax} v_{kl}
!! r_C^{n_{kl}} \exp(-dz_{kl} r_C)^2 \right) |l\rangle \langle l|
!!
double precision function Vps &
(a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl)
implicit none
integer n_a(3),n_b(3)
double precision g_a,g_b,a(3),b(3),c(3)
integer lmax,kmax,n_kl(kmax,0:lmax)
double precision v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
integer klocmax,n_k(klocmax)
double precision v_k(klocmax),dz_k(klocmax)
double precision Vloc,Vpseudo
Vps=Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) &
+Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c)
end
!!
!! Vps_num: brute force numerical evaluation of the same matrix element Vps
!!
double precision function Vps_num &
(npts,rmax,a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl)
implicit none
integer n_a(3),n_b(3)
double precision g_a,g_b,a(3),b(3),c(3),rmax
integer lmax,kmax,n_kl(kmax,0:lmax)
double precision v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
integer klocmax,n_k(klocmax)
double precision v_k(klocmax),dz_k(klocmax)
double precision Vloc_num,Vpseudo_num,v1,v2
integer npts,nptsgrid
nptsgrid=50
call initpseudos(nptsgrid)
v1=Vloc_num(npts,rmax,klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c)
v2=Vpseudo_num(nptsgrid,rmax,lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c)
Vps_num=v1+v2
end
double precision function Vloc_num(npts_over,xmax,klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c)
implicit none
integer klocmax
double precision v_k(klocmax),dz_k(klocmax)
integer n_k(klocmax)
integer npts_over,ix,iy,iz
double precision xmax,dx,x,y,z
double precision a(3),b(3),c(3),term,r,orb_phi,g_a,g_b,ac(3),bc(3)
integer n_a(3),n_b(3),k,l
do l=1,3
ac(l)=a(l)-c(l)
bc(l)=b(l)-c(l)
enddo
dx=2.d0*xmax/npts_over
Vloc_num=0.d0
do ix=1,npts_over
do iy=1,npts_over
do iz=1,npts_over
x=-xmax+dx*ix+dx/2.d0
y=-xmax+dx*iy+dx/2.d0
z=-xmax+dx*iz+dx/2.d0
term=orb_phi(x,y,z,n_a,ac,g_a)*orb_phi(x,y,z,n_b,bc,g_b)
r=dsqrt(x**2+y**2+z**2)
do k=1,klocmax
Vloc_num=Vloc_num+dx**3*v_k(k)*r**n_k(k)*dexp(-dz_k(k)*r**2)*term
enddo
enddo
enddo
enddo
end
double precision function orb_phi(x,y,z,npower,center,gamma)
implicit none
integer npower(3)
double precision x,y,z,r2,gamma,center(3)
r2=(x-center(1))**2+(y-center(2))**2+(z-center(3))**2
orb_phi=(x-center(1))**npower(1)*(y-center(2))**npower(2)*(z-center(3))**npower(3)
orb_phi=orb_phi*dexp(-gamma*r2)
end
!! Real spherical harmonics Ylm
! factor = ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2
! Y_lm(theta,phi) =
! m > 0 factor* P_l^|m|(cos(theta)) cos (|m| phi)
! m = 0 1/sqrt(2) *factor* P_l^0(cos(theta))
! m < 0 factor* P_l^|m|(cos(theta)) sin (|m| phi)
!
! x=cos(theta)
double precision function ylm_real(l,m,x,phi)
Using Intel IPP for sorting
2021-05-31 01:48:34 +02:00
implicit none
integer :: MM, iabs_m, m, l
double precision :: pi, fourpi, factor, x, phi, coef
double precision :: xchap, ychap, zchap
double precision, external :: fact
double precision :: PM(0:100,0:100), plm
Initial commit
2019-01-25 11:39:31 +01:00
MM=100
pi=dacos(-1.d0)
fourpi=4.d0*pi
iabs_m=iabs(m)
if(iabs_m.gt.l)stop 'm must be between -l and l'
factor= dsqrt( ((l+l+1)*fact(l-iabs_m))/(fourpi*fact(l+iabs_m)) )
if(dabs(x).gt.1.d0)then
print*,'pb. in ylm_no'
print*,'x=',x
stop
endif
call LPMN(MM,l,l,X,PM)
plm=PM(iabs_m,l)
coef=factor*plm
if(m.gt.0)ylm_real=dsqrt(2.d0)*coef*dcos(iabs_m*phi)
if(m.eq.0)ylm_real=coef
if(m.lt.0)ylm_real=dsqrt(2.d0)*coef*dsin(iabs_m*phi)
if(l.eq.0)ylm_real=dsqrt(1.d0/fourpi)
xchap=dsqrt(1.d0-x**2)*dcos(phi)
ychap=dsqrt(1.d0-x**2)*dsin(phi)
zchap=x
if(l.eq.1.and.m.eq.1)ylm_real=dsqrt(3.d0/fourpi)*xchap
if(l.eq.1.and.m.eq.0)ylm_real=dsqrt(3.d0/fourpi)*zchap
if(l.eq.1.and.m.eq.-1)ylm_real=dsqrt(3.d0/fourpi)*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.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
if(l.gt.2)stop 'l > 2 not coded!'
end
! _
! | |
! __ __ _ __ ___ ___ _ _ __| | ___
! \ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \
! \ V / | |_) \__ \ __/ |_| | (_| | (_) |
! \_/ | .__/|___/\___|\__,_|\____|\___/
! | |
! |_|
!! Routine Vpseudo is based on formumla (66)
!! of Kahn Baybutt TRuhlar J.Chem.Phys. vol.65 3826 (1976):
!!
!! Vpseudo= (4pi)**2* \sum_{l=0}^lmax \sum_{m=-l}^{l}
!! \sum{lambda=0}^{l+nA} \sum_{mu=-lambda}^{lambda}
!! \sum{k1=0}^{nAx} \sum{k2=0}^{nAy} \sum{k3=0}^{nAz}
!! binom(nAx,k1)*binom(nAy,k2)*binom(nAz,k3)* Y_{lambda mu}(AC_unit)
!! *CAx**(nAx-k1)*CAy**(nAy-k2)*CAz**(nAz-k3)*
!! bigI(lambda,mu,l,m,k1,k2,k3)
!! \sum{lambdap=0}^{l+nB} \sum_{mup=-lambdap}^{lambdap}
!! \sum{k1p=0}^{nBx} \sum{k2p=0}^{nBy} \sum{k3p=0}^{nBz}
!! binom(nBx,k1p)*binom(nBy,k2p)*binom(nBz,k3p)* Y_{lambdap mup}(BC_unit)
!! *CBx**(nBx-k1p)*CBy**(nBy-k2p)*CBz**(nBz-k3p)*
!! bigI(lambdap,mup,l,m,k1p,k2p,k3p)*
!! \sum_{k=1}^{kmax} v_kl(k,l)*
!! bigR(lambda,lambdap,k1+k2+k3+k1p+k2p+k3p+n_kl(k,l),g_a,g_b,AC,BC,dz_kl(k,l))
!!
!! nA=nAx+nAy+nAz
!! nB=nBx+nBy+nBz
!! AC=|A-C|
!! AC_x= A_x-C_x, etc.
!! BC=|B-C|
!! AC_unit= vect(AC)/AC
!! BC_unit= vect(BC)/BC
!! bigI(lambda,mu,l,m,k1,k2,k3)=
!! \int dOmega Y_{lambda mu}(Omega) xchap^k1 ychap^k2 zchap^k3 Y_{l m}(Omega)
!!
!! bigR(lambda,lambdap,N,g_a,g_b,gamm_k,AC,BC)
!! = exp(-g_a* AC**2 -g_b* BC**2) * int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal)
!! /int dx x^{ktot} exp(-g_k)*x**2) M_lambda(2 g_k D x)
double precision function Vpseudo &
(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c)
implicit none
! ___
! | ._ ._ _|_
! _|_ | | |_) |_| |_
! |
double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3)
integer, intent(in) :: lmax,kmax,n_kl(kmax,0:lmax)
integer, intent(in) :: n_a(3),n_b(3)
double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
!
! | _ _ _. |
! |_ (_) (_ (_| |
!
double precision :: fourpi,f,prod,prodp,binom_func,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
double precision :: arg
integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot
integer :: l,k, nkl_max
! _
! |_) o _ _. ._ ._ _.
! |_) | (_| (_| | | (_| \/
! _| /
double precision, allocatable :: array_coefs_A(:,:)
double precision, allocatable :: array_coefs_B(:,:)
double precision, allocatable :: array_R(:,:,:,:,:)
double precision, allocatable :: array_I_A(:,:,:,:,:)
double precision, allocatable :: array_I_B(:,:,:,:,:)
double precision :: f1, f2, f3
if (kmax.eq.1.and.lmax.eq.0.and.v_kl(1,0).eq.0.d0) then
Vpseudo=0.d0
return
end if
fourpi=4.d0*dacos(-1.d0)
ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2)
arg= g_a*ac*ac + g_b*bc*bc
if(arg.gt.-dlog(1.d-20))then
Vpseudo=0.d0
return
endif
freal=dexp(-arg)
areal=2.d0*g_a*ac
breal=2.d0*g_b*bc
ntotA=n_a(1)+n_a(2)+n_a(3)
ntotB=n_b(1)+n_b(2)+n_b(3)
ntot=ntotA+ntotB
nkl_max=4
allocate (array_coefs_A(0:ntot,3))
allocate (array_coefs_B(0:ntot,3))
allocate (array_R(kmax,0:ntot+nkl_max,0:lmax,0:lmax+ntot,0:lmax+ntot))
allocate (array_I_A(-(lmax+ntot):lmax+ntot,0:lmax+ntot,0:ntot,0:ntot,0:ntot))
allocate (array_I_B(-(lmax+ntot):lmax+ntot,0:lmax+ntot,0:ntot,0:ntot,0:ntot))
if(ac.eq.0.d0.and.bc.eq.0.d0)then
accu=0.d0
do k=1,kmax
do l=0,lmax
ktot=ntot+n_kl(k,l)
if (v_kl(k,l) == 0.d0) cycle
do m=-l,l
prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))
if (prod == 0.d0) cycle
prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))
if (prodp == 0.d0) cycle
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
Vpseudo=accu*fourpi
else if(ac.ne.0.d0.and.bc.ne.0.d0)then
f=fourpi*fourpi
theta_AC0=dacos( (a(3)-c(3))/ac )
phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac)
theta_BC0=dacos( (b(3)-c(3))/bc )
phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc)
do lambdap=0,lmax+ntotB
do lambda=0,lmax+ntotA
do l=0,lmax
do ktot=0,ntotA+ntotB+nkl_max
do k=1,kmax
array_R(k,ktot,l,lambda,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg)
enddo
enddo
enddo
enddo
enddo
do k1=0,n_a(1)
array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1)
enddo
do k2=0,n_a(2)
array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2)
enddo
do k3=0,n_a(3)
array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3)
enddo
do k1p=0,n_b(1)
array_coefs_B(k1p,1) = binom_func(n_b(1),k1p)*(c(1)-b(1))**(n_b(1)-k1p)
enddo
do k2p=0,n_b(2)
array_coefs_B(k2p,2) = binom_func(n_b(2),k2p)*(c(2)-b(2))**(n_b(2)-k2p)
enddo
do k3p=0,n_b(3)
array_coefs_B(k3p,3) = binom_func(n_b(3),k3p)*(c(3)-b(3))**(n_b(3)-k3p)
enddo
accu=0.d0
do l=0,lmax
do m=-l,l
do k3=0,n_a(3)
do k2=0,n_a(2)
do k1=0,n_a(1)
do lambda=0,l+ntotA
do mu=-lambda,lambda
array_I_A(mu,lambda,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
enddo
enddo
enddo
enddo
enddo
do k3p=0,n_b(3)
do k2p=0,n_b(2)
do k1p=0,n_b(1)
do lambdap=0,l+ntotB
do mup=-lambdap,lambdap
array_I_B(mup,lambdap,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p)
enddo
enddo
enddo
enddo
enddo
do k3=0,n_a(3)
if (array_coefs_A(k3,3) == 0.d0) cycle
do k2=0,n_a(2)
if (array_coefs_A(k2,2) == 0.d0) cycle
do k1=0,n_a(1)
if (array_coefs_A(k1,1) == 0.d0) cycle
do lambda=0,l+ntotA
do mu=-lambda,lambda
prod=ylm(lambda,mu,theta_AC0,phi_AC0)*array_coefs_A(k1,1)*array_coefs_A(k2,2)*array_coefs_A(k3,3)*array_I_A(mu,lambda,k1,k2,k3)
if (prod == 0.d0) cycle
do k3p=0,n_b(3)
do k2p=0,n_b(2)
do k1p=0,n_b(1)
do lambdap=0,l+ntotB
do mup=-lambdap,lambdap
prodp=prod*ylm(lambdap,mup,theta_BC0,phi_BC0)* &
array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)* &
array_I_B(mup,lambdap,k1p,k2p,k3p)
if (prodp == 0.d0) cycle
do k=1,kmax
ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l)
accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,lambdap)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
Vpseudo=f*accu
else if(ac.eq.0.d0.and.bc.ne.0.d0)then
f=fourpi**1.5d0
theta_BC0=dacos( (b(3)-c(3))/bc )
phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc)
areal=2.d0*g_a*ac
breal=2.d0*g_b*bc
freal=dexp(-g_a*ac**2-g_b*bc**2)
do lambdap=0,lmax+ntotB
do l=0,lmax
do ktot=0,ntotA+ntotB+nkl_max
do k=1,kmax
array_R(k,ktot,l,0,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
enddo
enddo
enddo
enddo
do k1p=0,n_b(1)
array_coefs_B(k1p,1) = binom_func(n_b(1),k1p)*(c(1)-b(1))**(n_b(1)-k1p)
enddo
do k2p=0,n_b(2)
array_coefs_B(k2p,2) = binom_func(n_b(2),k2p)*(c(2)-b(2))**(n_b(2)-k2p)
enddo
do k3p=0,n_b(3)
array_coefs_B(k3p,3) = binom_func(n_b(3),k3p)*(c(3)-b(3))**(n_b(3)-k3p)
enddo
accu=0.d0
do l=0,lmax
do m=-l,l
do k3p=0,n_b(3)
do k2p=0,n_b(2)
do k1p=0,n_b(1)
do lambdap=0,l+ntotB
do mup=-lambdap,lambdap
array_I_B(mup,lambdap,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p)
enddo
enddo
enddo
enddo
enddo
prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))
do k3p=0,n_b(3)
if (array_coefs_B(k3p,3) == 0.d0) cycle
do k2p=0,n_b(2)
if (array_coefs_B(k2p,2) == 0.d0) cycle
do k1p=0,n_b(1)
if (array_coefs_B(k1p,1) == 0.d0) cycle
do lambdap=0,l+ntotB
do mup=-lambdap,lambdap
prodp=prod*array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)*ylm(lambdap,mup,theta_BC0,phi_BC0)*array_I_B(mup,lambdap,k1p,k2p,k3p)
if (prodp == 0.d0) cycle
do k=1,kmax
ktot=ntotA+k1p+k2p+k3p+n_kl(k,l)
accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,0,lambdap)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
Vpseudo=f*accu
else if(ac.ne.0.d0.and.bc.eq.0.d0)then
f=fourpi**1.5d0
theta_AC0=dacos( (a(3)-c(3))/ac )
phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac)
areal=2.d0*g_a*ac
breal=2.d0*g_b*bc
freal=dexp(-g_a*ac**2-g_b*bc**2)
do lambda=0,lmax+ntotA
do l=0,lmax
do ktot=0,ntotA+ntotB+nkl_max
do k=1,kmax
array_R(k,ktot,l,lambda,0)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg)
enddo
enddo
enddo
enddo
do k1=0,n_a(1)
array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1)
enddo
do k2=0,n_a(2)
array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2)
enddo
do k3=0,n_a(3)
array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3)
enddo
accu=0.d0
do l=0,lmax
do m=-l,l
do k3=0,n_a(3)
do k2=0,n_a(2)
do k1=0,n_a(1)
do lambda=0,l+ntotA
do mu=-lambda,lambda
array_I_A(mu,lambda,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
enddo
enddo
enddo
enddo
enddo
do k3=0,n_a(3)
if (array_coefs_A(k3,3) == 0.d0) cycle
do k2=0,n_a(2)
if (array_coefs_A(k2,2) == 0.d0) cycle
do k1=0,n_a(1)
if (array_coefs_A(k1,1) == 0.d0) cycle
do lambda=0,l+ntotA
do mu=-lambda,lambda
prod=array_coefs_A(k1,1)*array_coefs_A(k2,2)*array_coefs_A(k3,3)*ylm(lambda,mu,theta_AC0,phi_AC0)*array_I_A(mu,lambda,k1,k2,k3)
if (prod == 0.d0) cycle
prodp=prod*bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))
if (prodp == 0.d0) cycle
do k=1,kmax
ktot=k1+k2+k3+ntotB+n_kl(k,l)
accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,0)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
Vpseudo=f*accu
endif
! _
! |_ o ._ _. | o _ _
! | | | | (_| | | _> (/_
!
deallocate (array_R, array_I_A, array_I_B)
deallocate (array_coefs_A, array_coefs_B)
return
end
! _
! | |
!__ __ _ __ ___ ___ _ _ __| | ___ _ __ _ _ _ __ ___
!\ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ | '_ \| | | | '_ ` _ \
! \ V / | |_) \__ \ __/ |_| | (_| | (_) | | | | | |_| | | | | | |
! \_/ | .__/|___/\___|\__,_|\__,_|\___/ |_| |_|\__,_|_| |_| |_|
! | |
! |_|
double precision function Vpseudo_num(npts,rmax,lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c)
implicit none
! ___
! | ._ ._ _|_
! _|_ | | |_) |_| |_
! |
double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3)
integer, intent(in) :: lmax,kmax,npts
integer, intent(in) :: n_a(3),n_b(3), n_kl(kmax,0:lmax)
double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
double precision, intent(in) :: rmax
!
! | _ _ _. |
! |_ (_) (_ (_| |
!
integer :: l,m,k,kk
double precision ac(3),bc(3)
double precision dr,sum,rC
double precision overlap_orb_ylm_brute
! _
! / _. | _ |
! \_ (_| | (_ |_| |
!
do l=1,3
ac(l)=a(l)-c(l)
bc(l)=b(l)-c(l)
enddo
dr=rmax/npts
sum=0.d0
do l=0,lmax
do m=-l,l
do k=1,npts
rC=(k-1)*dr+dr/2.d0
do kk=1,kmax
sum=sum+dr*v_kl(kk,l)*rC**(n_kl(kk,l)+2)*dexp(-dz_kl(kk,l)*rC**2) &
*overlap_orb_ylm_brute(npts,rC,n_a,ac,g_a,l,m) &
*overlap_orb_ylm_brute(npts,rC,n_b,bc,g_b,l,m)
enddo
enddo
enddo
enddo
Vpseudo_num=sum
return
end
!! Routine Vloc is a variation of formumla (66)
!! of Kahn Baybutt TRuhlar J.Chem.Phys. vol.65 3826 (1976)
!! without the projection operator
!!
!! Vloc= (4pi)**3/2* \sum_{k=1}^{klocmax} \sum_{l=0}^lmax \sum_{m=-l}^{l}
!!\sum{k1=0}^{nAx} \sum{k2=0}^{nAy} \sum{k3=0}^{nAz}
!! binom(nAx,k1)*binom(nAy,k2)*binom(nAz,k3)
!! *CAx**(nAx-k1)*CAy**(nAy-k2)*CAz**(nAz-k3)*
!! \sum{k1p=0}^{nBx} \sum{k2p=0}^{nBy} \sum{k3p=0}^{nBz}
!! binom(nBx,k1p)*binom(nBy,k2p)*binom(nBz,k3p)
!! *CBx**(nBx-k1p)*CBy**(nBy-k2p)*CBz**(nBz-k3p)*
!!\sum_{l=0}^lmax \sum_{m=-l}^{l}
!! bigI(0,0,l,m,k1+k1p,k2+k2p,k3+k3p)*Y_{l m}(D_unit)
!! *v_k(k)* bigR(lambda,k1+k2+k3+k1p+k2p+k3p+n_k(k),g_a,g_b,AC,BC,dz_k(k))
!!
!! nA=nAx+nAy+nAz
!! nB=nBx+nBy+nBz
!! D=(g_a AC+g_b BC)
!! D_unit= vect(D)/D
!! AC_x= A_x-C_x, etc.
!! BC=|B-C|
!! AC_unit= vect(AC)/AC
!! BC_unit= vect(BC)/BCA
!!
!! bigR(lambda,g_a,g_b,g_k,AC,BC)
!! = exp(-g_a* AC**2 -g_b* BC**2)*
!! I_loc= \int dx x**l *exp(-gam*x**2) M_n(ax) l=ktot+2 gam=g_a+g_b+dz_k(k) a=dreal n=l
!! M_n(x) modified spherical bessel function
double precision function Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c)
implicit none
integer klocmax
double precision v_k(klocmax),dz_k(klocmax),crochet,bigA
integer n_k(klocmax)
double precision a(3),g_a,b(3),g_b,c(3),d(3)
integer n_a(3),n_b(3),ntotA,ntotB,ntot,m
integer i,l,k,ktot,k1,k2,k3,k1p,k2p,k3p
double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0,coef
double precision,allocatable :: array_R_loc(:,:,:)
double precision,allocatable :: array_coefs(:,:,:,:,:,:)
double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg
fourpi=4.d0*dacos(-1.d0)
f=fourpi**1.5d0
ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2)
arg=g_a*ac**2+g_b*bc**2
Fix pseudos NaN
2021-07-06 15:02:47 +02:00
if(arg.gt.-dlog(1.d-20))then
Initial commit
2019-01-25 11:39:31 +01:00
Vloc=0.d0
return
endif
ntotA=n_a(1)+n_a(2)+n_a(3)
ntotB=n_b(1)+n_b(2)+n_b(3)
ntot=ntotA+ntotB
if(ac.eq.0.d0.and.bc.eq.0.d0)then
accu=0.d0
do k=1,klocmax
accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k))
enddo
Vloc=accu*fourpi*bigI(0,0,0,0,n_a(1)+n_b(1),n_a(2)+n_b(2),n_a(3)+n_b(3))
!bigI frequently is null
return
endif
freal=dexp(-g_a*ac**2-g_b*bc**2)
d2 = 0.d0
do i=1,3
d(i)=g_a*(a(i)-c(i))+g_b*(b(i)-c(i))
d2=d2+d(i)*d(i)
enddo
d2=dsqrt(d2)
dreal=2.d0*d2
allocate (array_R_loc(-2:ntot+klocmax,klocmax,0:ntot))
allocate (array_coefs(0:ntot,0:ntot,0:ntot,0:ntot,0:ntot,0:ntot))
do ktot=-2,ntotA+ntotB+klocmax
do l=0,ntot
do k=1,klocmax
array_R_loc(ktot,k,l)=freal*int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal)
enddo
enddo
enddo
do k1=0,n_a(1)
do k2=0,n_a(2)
do k3=0,n_a(3)
do k1p=0,n_b(1)
do k2p=0,n_b(2)
do k3p=0,n_b(3)
array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)&
*(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)&
*binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)&
*(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
enddo
enddo
enddo
enddo
enddo
enddo
accu=0.d0
if(d2 == 0.d0)then
l=0
m=0
coef=1.d0/dsqrt(4.d0*dacos(-1.d0))
do k=1,klocmax
do k1=0,n_a(1)
do k2=0,n_a(2)
do k3=0,n_a(3)
do k1p=0,n_b(1)
do k2p=0,n_b(2)
do k3p=0,n_b(3)
prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p) &
*bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p)
ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k)
accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
else
theta_DC0=dacos(d(3)/d2)
phi_DC0=datan2(d(2)/d2,d(1)/d2)
do k=1,klocmax
if (v_k(k) == 0.d0) cycle
do k1=0,n_a(1)
do k2=0,n_a(2)
do k3=0,n_a(3)
do k1p=0,n_b(1)
do k2p=0,n_b(2)
do k3p=0,n_b(3)
if (array_coefs(k1,k2,k3,k1p,k2p,k3p) == 0.d0) cycle
do l=0,ntot
do m=-l,l
coef=ylm(l,m,theta_DC0,phi_DC0)
if (coef == 0.d0) cycle
ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k)
if (array_R_loc(ktot,k,l) == 0.d0) cycle
prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p) &
*bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p)
accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
endif
Vloc=f*accu
deallocate (array_R_loc)
deallocate (array_coefs)
end
double precision function bigA(i,j,k)
implicit none
integer i,j,k
double precision fourpi,dble_fact
fourpi=4.d0*dacos(-1.d0)
bigA=0.d0
if(mod(i,2).eq.1)return
if(mod(j,2).eq.1)return
if(mod(k,2).eq.1)return
bigA=fourpi*dble_fact(i-1)*dble_fact(j-1)*dble_fact(k-1)/dble_fact(i+j+k+1)
end
!!
!! I_{lambda,mu,l,m}^{k1,k2,k3} = /int dOmega Y_{lambda mu} xchap^k1 ychap^k2 zchap^k3 Y_{lm}
!!
double precision function bigI(lambda,mu,l,m,k1,k2,k3)
implicit none
integer lambda,mu,l,m,k1,k2,k3
integer k,i,kp,ip
double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom_func,fact,coef_pm
double precision sgn, sgnp
pi=dacos(-1.d0)
bigI=0.d0
if(mu.gt.0.and.m.gt.0)then
sum=0.d0
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
if (factor1== 0.d0) return
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
if (factor2== 0.d0) return
sgn = 1.d0
do k=0,mu/2
do i=0,lambda-mu
if (coef_pm(lambda,i+mu) == 0.d0) cycle
sgnp = 1.d0
do kp=0,m/2
do ip=0,l-m
cylm=sgn*factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
if (cylm == 0.d0) cycle
cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
if (cylmp == 0.d0) cycle
sum=sum+cylm*cylmp*bigA(mu-2*k+m-2*kp+k1,2*k+2*kp+k2,i+ip+k3)
enddo
sgnp = -sgnp
enddo
enddo
sgn = -sgn
enddo
bigI=sum
return
endif
if(mu.eq.0.and.m.eq.0)then
factor1=dsqrt((2*lambda+1)/(4.d0*pi))
if (factor1== 0.d0) return
factor2=dsqrt((2*l+1)/(4.d0*pi))
if (factor2== 0.d0) return
sum=0.d0
do i=0,lambda
do ip=0,l
cylm=factor1*coef_pm(lambda,i)
if (cylm == 0.d0) cycle
cylmp=factor2*coef_pm(l,ip)
if (cylmp == 0.d0) cycle
sum=sum+cylm*cylmp*bigA(k1,k2,i+ip+k3)
enddo
enddo
bigI=sum
return
endif
if(mu.eq.0.and.m.gt.0)then
factor1=dsqrt((2*lambda+1)/(4.d0*pi))
if (factor1== 0.d0) return
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
if (factor2== 0.d0) return
sum=0.d0
do i=0,lambda
sgnp = 1.d0
do kp=0,m/2
do ip=0,l-m
cylm=factor1*coef_pm(lambda,i)
if (cylm == 0.d0) cycle
cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
if (cylmp == 0.d0) cycle
sum=sum+cylm*cylmp*bigA(m-2*kp+k1,2*kp+k2,i+ip+k3)
enddo
sgnp = -sgnp
enddo
enddo
bigI=sum
return
endif
if(mu.gt.0.and.m.eq.0)then
sum=0.d0
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
if (factor1== 0.d0) return
factor2=dsqrt((2*l+1)/(4.d0*pi))
if (factor2== 0.d0) return
sgn = 1.d0
do k=0,mu/2
do i=0,lambda-mu
if (coef_pm(lambda,i+mu) == 0.d0) cycle
do ip=0,l
cylm=sgn*factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
if (cylm == 0.d0) cycle
cylmp=factor2*coef_pm(l,ip)
if (cylmp == 0.d0) cycle
sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3)
enddo
enddo
sgn = -sgn
enddo
bigI=sum
return
endif
if(mu.lt.0.and.m.lt.0)then
mu=-mu
m=-m
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
if (factor1== 0.d0) return
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
if (factor2== 0.d0) return
sum=0.d0
sgn = 1.d0
do k=0,(mu-1)/2
do i=0,lambda-mu
if (coef_pm(lambda,i+mu) == 0.d0) cycle
sgnp = 1.d0
do kp=0,(m-1)/2
do ip=0,l-m
if (coef_pm(l,ip+m) == 0.d0) cycle
cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
if (cylm == 0.d0) cycle
cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
if (cylmp == 0.d0) cycle
sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-(2*kp+1)+k1,(2*k+1)+(2*kp+1)+k2,i+ip+k3)
enddo
sgnp = -sgnp
enddo
enddo
sgn = -sgn
enddo
mu=-mu
m=-m
bigI=sum
return
endif
if(mu.eq.0.and.m.lt.0)then
m=-m
factor1=dsqrt((2*lambda+1)/(4.d0*pi))
if (factor1 == 0.d0) return
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
if (factor2 == 0.d0) return
sum=0.d0
do i=0,lambda
sgnp = 1.d0
do kp=0,(m-1)/2
do ip=0,l-m
cylm=factor1*coef_pm(lambda,i)
if (cylm == 0.d0) cycle
cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
if (cylmp == 0.d0) cycle
sum=sum+cylm*cylmp*bigA(m-(2*kp+1)+k1,2*kp+1+k2,i+ip+k3)
enddo
sgnp = -sgnp
enddo
enddo
m=-m
bigI=sum
return
endif
if(mu.lt.0.and.m.eq.0)then
sum=0.d0
mu=-mu
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
if (factor1== 0.d0) return
factor2=dsqrt((2*l+1)/(4.d0*pi))
if (factor2== 0.d0) return
sgn = 1.d0
do k=0,(mu-1)/2
do i=0,lambda-mu
do ip=0,l
cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
if (cylm == 0.d0) cycle
cylmp=factor2*coef_pm(l,ip)
if (cylmp == 0.d0) cycle
sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3)
enddo
enddo
sgn = -sgn
enddo
mu=-mu
bigI=sum
return
endif
if(mu.gt.0.and.m.lt.0)then
sum=0.d0
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
if (factor1== 0.d0) return
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
if (factor2== 0.d0) return
m=-m
sgn=1.d0
do k=0,mu/2
do i=0,lambda-mu
if (coef_pm(lambda,i+mu) == 0.d0) cycle
sgnp=1.d0
do kp=0,(m-1)/2
do ip=0,l-m
if (coef_pm(l,ip+m) == 0.d0) cycle
cylm =sgn *factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
if (cylm == 0.d0) cycle
cylmp=sgnp*factor2*binom_func(m,