mirror of
https://github.com/LCPQ/quantum_package
synced 2024-09-13 06:38:31 +02:00
Accelerated pseudopotentials
This commit is contained in:
parent
1363b436e3
commit
5409607372
@ -66,6 +66,10 @@ Documentation
|
|||||||
title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
|
title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
|
||||||
|
|
||||||
|
|
||||||
|
`e_curve <http://github.com/LCPQ/quantum_package/tree/master/plugins/QmcChem/e_curve_qmc.irp.f#L1>`_
|
||||||
|
Undocumented
|
||||||
|
|
||||||
|
|
||||||
`mo_pseudo_grid <http://github.com/LCPQ/quantum_package/tree/master/plugins/QmcChem/pot_ao_pseudo_ints.irp.f#L56>`_
|
`mo_pseudo_grid <http://github.com/LCPQ/quantum_package/tree/master/plugins/QmcChem/pot_ao_pseudo_ints.irp.f#L56>`_
|
||||||
Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C
|
Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C
|
||||||
.br
|
.br
|
||||||
|
@ -350,6 +350,34 @@ subroutine write_spindeterminants
|
|||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, det_alpha_norm, (N_det_alpha_unique) ]
|
||||||
|
&BEGIN_PROVIDER [ double precision, det_beta_norm, (N_det_beta_unique) ]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Norm of the alpha and beta spin determinants in the wave function:
|
||||||
|
!
|
||||||
|
! ||Da||_i \sum_j C_{ij}**2
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
integer :: i,j,k,l
|
||||||
|
double precision :: f
|
||||||
|
|
||||||
|
det_alpha_norm = 0.d0
|
||||||
|
det_beta_norm = 0.d0
|
||||||
|
do k=1,N_det
|
||||||
|
i = psi_bilinear_matrix_rows(k)
|
||||||
|
j = psi_bilinear_matrix_columns(k)
|
||||||
|
do l=1,N_states
|
||||||
|
f = psi_bilinear_matrix_values(k,l)*psi_bilinear_matrix_values(k,l)
|
||||||
|
enddo
|
||||||
|
det_alpha_norm(i) += f
|
||||||
|
det_beta_norm(j) += f
|
||||||
|
enddo
|
||||||
|
det_alpha_norm = det_alpha_norm / dble(N_states)
|
||||||
|
det_beta_norm = det_beta_norm / dble(N_states)
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
!==============================================================================!
|
!==============================================================================!
|
||||||
! !
|
! !
|
||||||
|
@ -214,13 +214,14 @@ integer :: l,k, nkl_max
|
|||||||
! |_) | (_| (_| | | (_| \/
|
! |_) | (_| (_| | | (_| \/
|
||||||
! _| /
|
! _| /
|
||||||
|
|
||||||
double precision, allocatable :: array_coefs_A(:,:,:)
|
double precision, allocatable :: array_coefs_A(:,:)
|
||||||
double precision, allocatable :: array_coefs_B(:,:,:)
|
double precision, allocatable :: array_coefs_B(:,:)
|
||||||
|
|
||||||
double precision, allocatable :: array_R(:,:,:,:,:)
|
double precision, allocatable :: array_R(:,:,:,:,:)
|
||||||
double precision, allocatable :: array_I_A(:,:,:,:,:)
|
double precision, allocatable :: array_I_A(:,:,:,:,:)
|
||||||
double precision, allocatable :: array_I_B(:,:,:,:,:)
|
double precision, allocatable :: array_I_B(:,:,:,:,:)
|
||||||
|
|
||||||
|
double precision :: f1, f2, f3
|
||||||
|
|
||||||
! _
|
! _
|
||||||
! / _. | _ |
|
! / _. | _ |
|
||||||
@ -255,14 +256,14 @@ nkl_max=4
|
|||||||
! A l l o c a t e !
|
! A l l o c a t e !
|
||||||
!=!=!=!=!=!=!=!=!=!
|
!=!=!=!=!=!=!=!=!=!
|
||||||
|
|
||||||
allocate (array_coefs_A(0:ntot,0:ntot,0:ntot))
|
allocate (array_coefs_A(0:ntot,3))
|
||||||
allocate (array_coefs_B(0:ntot,0:ntot,0:ntot))
|
allocate (array_coefs_B(0:ntot,3))
|
||||||
|
|
||||||
allocate (array_R(0:ntot+nkl_max,kmax,0:lmax,0:lmax+ntot,0:lmax+ntot))
|
allocate (array_R(kmax,0:ntot+nkl_max,0:lmax,0:lmax+ntot,0:lmax+ntot))
|
||||||
|
|
||||||
allocate (array_I_A(0:lmax+ntot,-(lmax+ntot):lmax+ntot,0:ntot,0:ntot,0:ntot))
|
allocate (array_I_A(-(lmax+ntot):lmax+ntot,0:lmax+ntot,0:ntot,0:ntot,0:ntot))
|
||||||
|
|
||||||
allocate (array_I_B(0:lmax+ntot,-(lmax+ntot):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
|
if(ac.eq.0.d0.and.bc.eq.0.d0)then
|
||||||
|
|
||||||
@ -310,14 +311,12 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then
|
|||||||
phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc)
|
phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
do ktot=0,ntotA+ntotB+nkl_max
|
|
||||||
do lambda=0,lmax+ntotA
|
|
||||||
do lambdap=0,lmax+ntotB
|
do lambdap=0,lmax+ntotB
|
||||||
do k=1,kmax
|
do lambda=0,lmax+ntotA
|
||||||
do l=0,lmax
|
do l=0,lmax
|
||||||
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)
|
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
|
enddo
|
||||||
@ -325,21 +324,23 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then
|
|||||||
enddo
|
enddo
|
||||||
|
|
||||||
do k1=0,n_a(1)
|
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)
|
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)
|
do k3=0,n_a(3)
|
||||||
array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3) &
|
array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(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)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do k1p=0,n_b(1)
|
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)
|
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)
|
do k3p=0,n_b(3)
|
||||||
array_coefs_B(k1p,k2p,k3p)=binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p) &
|
array_coefs_B(k3p,3) = binom_func(n_b(3),k3p)*(c(3)-b(3))**(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
|
||||||
|
|
||||||
!=!=!=!=!=!=!=!
|
!=!=!=!=!=!=!=!
|
||||||
@ -350,51 +351,53 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then
|
|||||||
do l=0,lmax
|
do l=0,lmax
|
||||||
do m=-l,l
|
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 lambda=0,l+ntotA
|
||||||
do mu=-lambda,lambda
|
do mu=-lambda,lambda
|
||||||
do k1=0,n_a(1)
|
array_I_A(mu,lambda,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
|
||||||
do k2=0,n_a(2)
|
|
||||||
do k3=0,n_a(3)
|
|
||||||
array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
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 lambdap=0,l+ntotB
|
||||||
do mup=-lambdap,lambdap
|
do mup=-lambdap,lambdap
|
||||||
do k1p=0,n_b(1)
|
array_I_B(mup,lambdap,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p)
|
||||||
do k2p=0,n_b(2)
|
|
||||||
do k3p=0,n_b(3)
|
|
||||||
array_I_B(lambdap,mup,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p)
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
|
do k3=0,n_a(3)
|
||||||
|
do k2=0,n_a(2)
|
||||||
|
do k1=0,n_a(1)
|
||||||
|
|
||||||
do lambda=0,l+ntotA
|
do lambda=0,l+ntotA
|
||||||
do mu=-lambda,lambda
|
do mu=-lambda,lambda
|
||||||
|
|
||||||
do k1=0,n_a(1)
|
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)
|
||||||
do k2=0,n_a(2)
|
|
||||||
do k3=0,n_a(3)
|
|
||||||
|
|
||||||
prod=ylm(lambda,mu,theta_AC0,phi_AC0)*array_coefs_A(k1,k2,k3)*array_I_A(lambda,mu,k1,k2,k3)
|
|
||||||
|
|
||||||
|
do k3p=0,n_b(3)
|
||||||
|
do k2p=0,n_b(2)
|
||||||
|
do k1p=0,n_b(1)
|
||||||
do lambdap=0,l+ntotB
|
do lambdap=0,l+ntotB
|
||||||
do mup=-lambdap,lambdap
|
do mup=-lambdap,lambdap
|
||||||
|
|
||||||
do k1p=0,n_b(1)
|
prodp=prod*ylm(lambdap,mup,theta_BC0,phi_BC0)* &
|
||||||
do k2p=0,n_b(2)
|
array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)* &
|
||||||
do k3p=0,n_b(3)
|
array_I_B(mup,lambdap,k1p,k2p,k3p)
|
||||||
|
|
||||||
prodp=ylm(lambdap,mup,theta_BC0,phi_BC0)*array_coefs_B(k1p,k2p,k3p)*array_I_B(lambdap,mup,k1p,k2p,k3p)
|
|
||||||
|
|
||||||
do k=1,kmax
|
do k=1,kmax
|
||||||
ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l)
|
ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l)
|
||||||
accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,lambda,lambdap)
|
accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,lambdap)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
@ -434,24 +437,24 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
|
|||||||
breal=2.d0*g_b*bc
|
breal=2.d0*g_b*bc
|
||||||
freal=dexp(-g_a*ac**2-g_b*bc**2)
|
freal=dexp(-g_a*ac**2-g_b*bc**2)
|
||||||
|
|
||||||
do ktot=0,ntotA+ntotB+nkl_max
|
|
||||||
do lambdap=0,lmax+ntotB
|
do lambdap=0,lmax+ntotB
|
||||||
do k=1,kmax
|
|
||||||
do l=0,lmax
|
do l=0,lmax
|
||||||
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)
|
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
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do k1p=0,n_b(1)
|
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)
|
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)
|
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)
|
||||||
array_coefs_B(k1p,k2p,k3p)=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
|
||||||
|
|
||||||
!=!=!=!=!=!=!=!
|
!=!=!=!=!=!=!=!
|
||||||
@ -462,12 +465,12 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
|
|||||||
do l=0,lmax
|
do l=0,lmax
|
||||||
do m=-l,l
|
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 lambdap=0,l+ntotB
|
||||||
do mup=-lambdap,lambdap
|
do mup=-lambdap,lambdap
|
||||||
do k1p=0,n_b(1)
|
array_I_B(mup,lambdap,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p)
|
||||||
do k2p=0,n_b(2)
|
|
||||||
do k3p=0,n_b(3)
|
|
||||||
array_I_B(lambdap,mup,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p)
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -476,18 +479,18 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
|
|||||||
|
|
||||||
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))
|
||||||
|
|
||||||
|
do k3p=0,n_b(3)
|
||||||
|
do k2p=0,n_b(2)
|
||||||
|
do k1p=0,n_b(1)
|
||||||
do lambdap=0,l+ntotB
|
do lambdap=0,l+ntotB
|
||||||
do mup=-lambdap,lambdap
|
do mup=-lambdap,lambdap
|
||||||
do k1p=0,n_b(1)
|
|
||||||
do k2p=0,n_b(2)
|
|
||||||
do k3p=0,n_b(3)
|
|
||||||
|
|
||||||
prodp=array_coefs_B(k1p,k2p,k3p)*ylm(lambdap,mup,theta_BC0,phi_BC0)*array_I_B(lambdap,mup,k1p,k2p,k3p)
|
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)
|
||||||
|
|
||||||
do k=1,kmax
|
do k=1,kmax
|
||||||
|
|
||||||
ktot=ntotA+k1p+k2p+k3p+n_kl(k,l)
|
ktot=ntotA+k1p+k2p+k3p+n_kl(k,l)
|
||||||
accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,0,lambdap)
|
accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,0,lambdap)
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
@ -519,26 +522,24 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then
|
|||||||
breal=2.d0*g_b*bc
|
breal=2.d0*g_b*bc
|
||||||
freal=dexp(-g_a*ac**2-g_b*bc**2)
|
freal=dexp(-g_a*ac**2-g_b*bc**2)
|
||||||
|
|
||||||
do ktot=0,ntotA+ntotB+nkl_max
|
|
||||||
do lambda=0,lmax+ntotA
|
do lambda=0,lmax+ntotA
|
||||||
do k=1,kmax
|
|
||||||
do l=0,lmax
|
do l=0,lmax
|
||||||
|
do ktot=0,ntotA+ntotB+nkl_max
|
||||||
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)
|
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
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do k1=0,n_a(1)
|
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)
|
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)
|
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)
|
||||||
array_coefs_A(k1,k2,k3)=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)
|
|
||||||
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
!=!=!=!=!=!=!=!
|
!=!=!=!=!=!=!=!
|
||||||
@ -549,30 +550,30 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then
|
|||||||
do l=0,lmax
|
do l=0,lmax
|
||||||
do m=-l,l
|
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 lambda=0,l+ntotA
|
||||||
do mu=-lambda,lambda
|
do mu=-lambda,lambda
|
||||||
do k1=0,n_a(1)
|
array_I_A(mu,lambda,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
|
||||||
do k2=0,n_a(2)
|
|
||||||
do k3=0,n_a(3)
|
|
||||||
array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
|
do k3=0,n_a(3)
|
||||||
|
do k2=0,n_a(2)
|
||||||
|
do k1=0,n_a(1)
|
||||||
do lambda=0,l+ntotA
|
do lambda=0,l+ntotA
|
||||||
do mu=-lambda,lambda
|
do mu=-lambda,lambda
|
||||||
do k1=0,n_a(1)
|
|
||||||
do k2=0,n_a(2)
|
|
||||||
do k3=0,n_a(3)
|
|
||||||
|
|
||||||
prod=array_coefs_A(k1,k2,k3)*ylm(lambda,mu,theta_AC0,phi_AC0)*array_I_A(lambda,mu,k1,k2,k3)
|
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)
|
||||||
prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))
|
prodp=prod*bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))
|
||||||
|
|
||||||
do k=1,kmax
|
do k=1,kmax
|
||||||
ktot=k1+k2+k3+ntotB+n_kl(k,l)
|
ktot=k1+k2+k3+ntotB+n_kl(k,l)
|
||||||
accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,lambda,0)
|
accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,0)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
@ -850,22 +851,27 @@ implicit none
|
|||||||
integer lambda,mu,l,m,k1,k2,k3
|
integer lambda,mu,l,m,k1,k2,k3
|
||||||
integer k,i,kp,ip
|
integer k,i,kp,ip
|
||||||
double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom_func,fact,coef_pm
|
double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom_func,fact,coef_pm
|
||||||
|
double precision sgn, sgnp
|
||||||
pi=dacos(-1.d0)
|
pi=dacos(-1.d0)
|
||||||
|
|
||||||
if(mu.gt.0.and.m.gt.0)then
|
if(mu.gt.0.and.m.gt.0)then
|
||||||
sum=0.d0
|
sum=0.d0
|
||||||
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
|
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
|
||||||
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
|
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
|
||||||
|
sgn = 1.d0
|
||||||
do k=0,mu/2
|
do k=0,mu/2
|
||||||
do i=0,lambda-mu
|
do i=0,lambda-mu
|
||||||
|
sgnp = 1.d0
|
||||||
do kp=0,m/2
|
do kp=0,m/2
|
||||||
do ip=0,l-m
|
do ip=0,l-m
|
||||||
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
cylm=sgn*factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
||||||
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
||||||
sum=sum+cylm*cylmp*bigA(mu-2*k+m-2*kp+k1,2*k+2*kp+k2,i+ip+k3)
|
sum=sum+cylm*cylmp*bigA(mu-2*k+m-2*kp+k1,2*k+2*kp+k2,i+ip+k3)
|
||||||
enddo
|
enddo
|
||||||
|
sgnp = -sgnp
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
sgn = -sgn
|
||||||
enddo
|
enddo
|
||||||
bigI=sum
|
bigI=sum
|
||||||
return
|
return
|
||||||
@ -888,15 +894,17 @@ endif
|
|||||||
|
|
||||||
if(mu.eq.0.and.m.gt.0)then
|
if(mu.eq.0.and.m.gt.0)then
|
||||||
factor1=dsqrt((2*lambda+1)/(4.d0*pi))
|
factor1=dsqrt((2*lambda+1)/(4.d0*pi))
|
||||||
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
|
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
|
||||||
sum=0.d0
|
sum=0.d0
|
||||||
do i=0,lambda
|
do i=0,lambda
|
||||||
|
sgnp = 1.d0
|
||||||
do kp=0,m/2
|
do kp=0,m/2
|
||||||
do ip=0,l-m
|
do ip=0,l-m
|
||||||
cylm=factor1*coef_pm(lambda,i)
|
cylm=factor1*coef_pm(lambda,i)
|
||||||
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
||||||
sum=sum+cylm*cylmp*bigA(m-2*kp+k1,2*kp+k2,i+ip+k3)
|
sum=sum+cylm*cylmp*bigA(m-2*kp+k1,2*kp+k2,i+ip+k3)
|
||||||
enddo
|
enddo
|
||||||
|
sgnp = -sgnp
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
bigI=sum
|
bigI=sum
|
||||||
@ -905,16 +913,18 @@ endif
|
|||||||
|
|
||||||
if(mu.gt.0.and.m.eq.0)then
|
if(mu.gt.0.and.m.eq.0)then
|
||||||
sum=0.d0
|
sum=0.d0
|
||||||
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
|
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
|
||||||
factor2=dsqrt((2*l+1)/(4.d0*pi))
|
factor2=dsqrt((2*l+1)/(4.d0*pi))
|
||||||
|
sgn = 1.d0
|
||||||
do k=0,mu/2
|
do k=0,mu/2
|
||||||
do i=0,lambda-mu
|
do i=0,lambda-mu
|
||||||
do ip=0,l
|
do ip=0,l
|
||||||
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
cylm=sgn*factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
||||||
cylmp=factor2*coef_pm(l,ip)
|
cylmp=factor2*coef_pm(l,ip)
|
||||||
sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3)
|
sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
sgn = -sgn
|
||||||
enddo
|
enddo
|
||||||
bigI=sum
|
bigI=sum
|
||||||
return
|
return
|
||||||
@ -923,19 +933,23 @@ endif
|
|||||||
if(mu.lt.0.and.m.lt.0)then
|
if(mu.lt.0.and.m.lt.0)then
|
||||||
mu=-mu
|
mu=-mu
|
||||||
m=-m
|
m=-m
|
||||||
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
|
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
|
||||||
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
|
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
|
||||||
sum=0.d0
|
sum=0.d0
|
||||||
|
sgn = 1.d0
|
||||||
do k=0,(mu-1)/2
|
do k=0,(mu-1)/2
|
||||||
do i=0,lambda-mu
|
do i=0,lambda-mu
|
||||||
|
sgnp = 1.d0
|
||||||
do kp=0,(m-1)/2
|
do kp=0,(m-1)/2
|
||||||
do ip=0,l-m
|
do ip=0,l-m
|
||||||
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
||||||
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
||||||
sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-(2*kp+1)+k1,(2*k+1)+(2*kp+1)+k2,i+ip+k3)
|
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
|
enddo
|
||||||
|
sgnp = -sgnp
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
sgn = -sgn
|
||||||
enddo
|
enddo
|
||||||
mu=-mu
|
mu=-mu
|
||||||
m=-m
|
m=-m
|
||||||
@ -946,15 +960,17 @@ endif
|
|||||||
if(mu.eq.0.and.m.lt.0)then
|
if(mu.eq.0.and.m.lt.0)then
|
||||||
m=-m
|
m=-m
|
||||||
factor1=dsqrt((2*lambda+1)/(4.d0*pi))
|
factor1=dsqrt((2*lambda+1)/(4.d0*pi))
|
||||||
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
|
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
|
||||||
sum=0.d0
|
sum=0.d0
|
||||||
do i=0,lambda
|
do i=0,lambda
|
||||||
|
sgnp = 1.d0
|
||||||
do kp=0,(m-1)/2
|
do kp=0,(m-1)/2
|
||||||
do ip=0,l-m
|
do ip=0,l-m
|
||||||
cylm=factor1*coef_pm(lambda,i)
|
cylm=factor1*coef_pm(lambda,i)
|
||||||
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
||||||
sum=sum+cylm*cylmp*bigA(m-(2*kp+1)+k1,2*kp+1+k2,i+ip+k3)
|
sum=sum+cylm*cylmp*bigA(m-(2*kp+1)+k1,2*kp+1+k2,i+ip+k3)
|
||||||
enddo
|
enddo
|
||||||
|
sgnp = -sgnp
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
m=-m
|
m=-m
|
||||||
@ -965,16 +981,18 @@ endif
|
|||||||
if(mu.lt.0.and.m.eq.0)then
|
if(mu.lt.0.and.m.eq.0)then
|
||||||
sum=0.d0
|
sum=0.d0
|
||||||
mu=-mu
|
mu=-mu
|
||||||
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
|
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
|
||||||
factor2=dsqrt((2*l+1)/(4.d0*pi))
|
factor2=dsqrt((2*l+1)/(4.d0*pi))
|
||||||
|
sgn = 1.d0
|
||||||
do k=0,(mu-1)/2
|
do k=0,(mu-1)/2
|
||||||
do i=0,lambda-mu
|
do i=0,lambda-mu
|
||||||
do ip=0,l
|
do ip=0,l
|
||||||
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
||||||
cylmp=factor2*coef_pm(l,ip)
|
cylmp=factor2*coef_pm(l,ip)
|
||||||
sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3)
|
sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
sgn = -sgn
|
||||||
enddo
|
enddo
|
||||||
mu=-mu
|
mu=-mu
|
||||||
bigI=sum
|
bigI=sum
|
||||||
@ -983,19 +1001,23 @@ endif
|
|||||||
|
|
||||||
if(mu.gt.0.and.m.lt.0)then
|
if(mu.gt.0.and.m.lt.0)then
|
||||||
sum=0.d0
|
sum=0.d0
|
||||||
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
|
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
|
||||||
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
|
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
|
||||||
m=-m
|
m=-m
|
||||||
|
sgn=1.d0
|
||||||
do k=0,mu/2
|
do k=0,mu/2
|
||||||
do i=0,lambda-mu
|
do i=0,lambda-mu
|
||||||
|
sgnp=1.d0
|
||||||
do kp=0,(m-1)/2
|
do kp=0,(m-1)/2
|
||||||
do ip=0,l-m
|
do ip=0,l-m
|
||||||
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
cylm =sgn *factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
||||||
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
||||||
sum=sum+cylm*cylmp*bigA(mu-2*k+m-(2*kp+1)+k1,2*k+2*kp+1+k2,i+ip+k3)
|
sum=sum+cylm*cylmp*bigA(mu-2*k+m-(2*kp+1)+k1,2*k+2*kp+1+k2,i+ip+k3)
|
||||||
enddo
|
enddo
|
||||||
|
sgnp = -sgnp
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
sgn = -sgn
|
||||||
enddo
|
enddo
|
||||||
m=-m
|
m=-m
|
||||||
bigI=sum
|
bigI=sum
|
||||||
@ -1004,19 +1026,23 @@ endif
|
|||||||
|
|
||||||
if(mu.lt.0.and.m.gt.0)then
|
if(mu.lt.0.and.m.gt.0)then
|
||||||
mu=-mu
|
mu=-mu
|
||||||
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
|
factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
|
||||||
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
|
factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
|
||||||
sum=0.d0
|
sum=0.d0
|
||||||
|
sgn = 1.d0
|
||||||
do k=0,(mu-1)/2
|
do k=0,(mu-1)/2
|
||||||
do i=0,lambda-mu
|
do i=0,lambda-mu
|
||||||
|
sgnp = 1.d0
|
||||||
do kp=0,m/2
|
do kp=0,m/2
|
||||||
do ip=0,l-m
|
do ip=0,l-m
|
||||||
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
cylm=sgn*factor1 *binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
|
||||||
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
|
||||||
sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-2*kp+k1,2*k+1+2*kp+k2,i+ip+k3)
|
sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-2*kp+k1,2*k+1+2*kp+k2,i+ip+k3)
|
||||||
enddo
|
enddo
|
||||||
|
sgnp = -sgnp
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
sgn = -sgn
|
||||||
enddo
|
enddo
|
||||||
bigI=sum
|
bigI=sum
|
||||||
mu=-mu
|
mu=-mu
|
||||||
@ -1128,28 +1154,49 @@ end
|
|||||||
!
|
!
|
||||||
IMPLICIT DOUBLE PRECISION (P,X)
|
IMPLICIT DOUBLE PRECISION (P,X)
|
||||||
DIMENSION PM(0:MM,0:(N+1))
|
DIMENSION PM(0:MM,0:(N+1))
|
||||||
DO 10 I=0,N
|
DOUBLE PRECISION, SAVE :: INVERSE(100) = 0.D0
|
||||||
DO 10 J=0,M
|
DOUBLE PRECISION :: LS, II, JJ
|
||||||
10 PM(J,I)=0.0D0
|
IF (INVERSE(1) == 0.d0) THEN
|
||||||
|
DO I=1,100
|
||||||
|
INVERSE(I) = 1.D0/DBLE(I)
|
||||||
|
ENDDO
|
||||||
|
ENDIF
|
||||||
|
DO I=0,N
|
||||||
|
DO J=0,M
|
||||||
|
PM(J,I)=0.0D0
|
||||||
|
ENDDO
|
||||||
|
ENDDO
|
||||||
PM(0,0)=1.0D0
|
PM(0,0)=1.0D0
|
||||||
IF (DABS(X).EQ.1.0D0) THEN
|
IF (DABS(X).EQ.1.0D0) THEN
|
||||||
DO 15 I=1,N
|
DO I=1,N
|
||||||
15 PM(0,I)=X**I
|
PM(0,I)=X**I
|
||||||
|
ENDDO
|
||||||
RETURN
|
RETURN
|
||||||
ENDIF
|
ENDIF
|
||||||
LS=1
|
LS=1.D0
|
||||||
IF (DABS(X).GT.1.0D0) LS=-1
|
IF (DABS(X).GT.1.0D0) LS=-1.D0
|
||||||
XQ=DSQRT(LS*(1.0D0-X*X))
|
XQ=DSQRT(LS*(1.0D0-X*X))
|
||||||
XS=LS*(1.0D0-X*X)
|
XS=LS*(1.0D0-X*X)
|
||||||
DO 30 I=1,M
|
II = 1.D0
|
||||||
30 PM(I,I)=-LS*(2.0D0*I-1.0D0)*XQ*PM(I-1,I-1)
|
DO I=1,M
|
||||||
DO 35 I=0,M
|
PM(I,I)=-LS*II*XQ*PM(I-1,I-1)
|
||||||
35 PM(I,I+1)=(2.0D0*I+1.0D0)*X*PM(I,I)
|
II = II+2.D0
|
||||||
|
ENDDO
|
||||||
|
II = 1.D0
|
||||||
|
DO I=0,M
|
||||||
|
PM(I,I+1)=II*X*PM(I,I)
|
||||||
|
II = II+2.D0
|
||||||
|
ENDDO
|
||||||
|
|
||||||
DO 40 I=0,M
|
II = 0.D0
|
||||||
DO 40 J=I+2,N
|
DO I=0,M
|
||||||
PM(I,J)=((2.0D0*J-1.0D0)*X*PM(I,J-1)- (I+J-1.0D0)*PM(I,J-2))/(J-I)
|
JJ = II+2.D0
|
||||||
40 CONTINUE
|
DO J=I+2,N
|
||||||
|
PM(I,J)=((2.0D0*JJ-1.0D0)*X*PM(I,J-1)- (II+JJ-1.0D0)*PM(I,J-2))*INVERSE(J-I)
|
||||||
|
JJ = JJ+1.D0
|
||||||
|
ENDDO
|
||||||
|
II = II+1.D0
|
||||||
|
ENDDO
|
||||||
END
|
END
|
||||||
|
|
||||||
|
|
||||||
@ -1703,17 +1750,37 @@ end
|
|||||||
!c
|
!c
|
||||||
double precision function ylm(l,m,theta,phi)
|
double precision function ylm(l,m,theta,phi)
|
||||||
implicit none
|
implicit none
|
||||||
integer l,m
|
integer l,m,i
|
||||||
double precision theta,phi,pm,factor,pi,x,fact,sign
|
double precision theta,phi,pm,factor,twopi,x,fact,sign
|
||||||
DIMENSION PM(0:100,0:100)
|
DIMENSION PM(0:100,0:100)
|
||||||
pi=dacos(-1.d0)
|
twopi=2.d0*dacos(-1.d0)
|
||||||
x=dcos(theta)
|
x=dcos(theta)
|
||||||
sign=(-1.d0)**m
|
if (iand(m,1)) then
|
||||||
|
sign=-1.d0
|
||||||
|
else
|
||||||
|
sign=1.d0
|
||||||
|
endif
|
||||||
CALL LPMN(100,l,l,X,PM)
|
CALL LPMN(100,l,l,X,PM)
|
||||||
factor=dsqrt( (2*l+1)*fact(l-iabs(m)) /(4.d0*pi*fact(l+iabs(m))) )
|
if (m > 0) then
|
||||||
if(m.gt.0)ylm=sign*dsqrt(2.d0)*factor*pm(m,l)*dcos(dfloat(m)*phi)
|
factor=dsqrt((l+l+1)*fact(l-m) /(twopi*fact(l+m)) )
|
||||||
if(m.eq.0)ylm=factor*pm(m,l)
|
! factor = dble(l+m)
|
||||||
if(m.lt.0)ylm=sign*dsqrt(2.d0)*factor*pm(iabs(m),l)*dsin(dfloat(iabs(m))*phi)
|
! do i=-m,m-1
|
||||||
|
! factor = factor * (factor - 1.d0)
|
||||||
|
! enddo
|
||||||
|
! factor=dsqrt(dble(l+l+1)/(twopi*factor) )
|
||||||
|
ylm=sign*factor*pm(m,l)*dcos(dfloat(m)*phi)
|
||||||
|
else if (m == 0) then
|
||||||
|
factor=dsqrt( 0.5d0*(l+l+1) /twopi )
|
||||||
|
ylm=factor*pm(m,l)
|
||||||
|
else if (m < 0) then
|
||||||
|
factor=dsqrt( (l+l+1)*fact(l+m) /(twopi*fact(l-m)) )
|
||||||
|
! factor = dble(l-m)
|
||||||
|
! do i=m,-m-1
|
||||||
|
! factor = factor * (factor - 1.d0)
|
||||||
|
! enddo
|
||||||
|
! factor=dsqrt(dble(l+l+1)/(twopi*factor) )
|
||||||
|
ylm=sign*factor*pm(-m,l)*dsin(dfloat(-m)*phi)
|
||||||
|
endif
|
||||||
end
|
end
|
||||||
|
|
||||||
!c Explicit representation of Legendre polynomials P_n(x)
|
!c Explicit representation of Legendre polynomials P_n(x)
|
||||||
@ -2139,10 +2206,11 @@ double precision prod,binom_func,accu,bigI,ylm,bessel_mod
|
|||||||
double precision theta_AC0,phi_AC0,ac,ac2,factor,fourpi,arg,r,areal
|
double precision theta_AC0,phi_AC0,ac,ac2,factor,fourpi,arg,r,areal
|
||||||
integer ntotA,mu,k1,k2,k3,lambda
|
integer ntotA,mu,k1,k2,k3,lambda
|
||||||
integer n_a(3)
|
integer n_a(3)
|
||||||
double precision array_coefs_A(0:ntot_max,0:ntot_max,0:ntot_max), y
|
double precision y, f1, f2
|
||||||
|
double precision, allocatable :: array_coefs_A(:,:)
|
||||||
|
|
||||||
ac2=(a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2
|
ac2=(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)
|
ac=dsqrt(ac2)
|
||||||
arg=g_a*(ac2+r*r)
|
arg=g_a*(ac2+r*r)
|
||||||
fourpi=4.d0*dacos(-1.d0)
|
fourpi=4.d0*dacos(-1.d0)
|
||||||
factor=fourpi*dexp(-arg)
|
factor=fourpi*dexp(-arg)
|
||||||
@ -2159,14 +2227,15 @@ else
|
|||||||
theta_AC0=dacos( (a(3)-c(3))/ac )
|
theta_AC0=dacos( (a(3)-c(3))/ac )
|
||||||
phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac)
|
phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac)
|
||||||
|
|
||||||
|
allocate (array_coefs_A(0:ntotA,3))
|
||||||
do k1=0,n_a(1)
|
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)*r**(k1)
|
||||||
|
enddo
|
||||||
do k2=0,n_a(2)
|
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)*r**(k2)
|
||||||
|
enddo
|
||||||
do k3=0,n_a(3)
|
do k3=0,n_a(3)
|
||||||
array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3) &
|
array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3)*r**(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) &
|
|
||||||
*r**(k1+k2+k3)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
accu=0.d0
|
accu=0.d0
|
||||||
@ -2176,10 +2245,14 @@ else
|
|||||||
if (y == 0.d0) then
|
if (y == 0.d0) then
|
||||||
cycle
|
cycle
|
||||||
endif
|
endif
|
||||||
do k1=0,n_a(1)
|
|
||||||
do k2=0,n_a(2)
|
|
||||||
do k3=0,n_a(3)
|
do k3=0,n_a(3)
|
||||||
prod=y*array_coefs_A(k1,k2,k3)*bigI(lambda,mu,l,m,k1,k2,k3)
|
f1 = y*array_coefs_A(k3,3)
|
||||||
|
if (f1 == 0.d0) cycle
|
||||||
|
do k2=0,n_a(2)
|
||||||
|
f2 = f1*array_coefs_A(k2,2)
|
||||||
|
if (f2 == 0.d0) cycle
|
||||||
|
do k1=0,n_a(1)
|
||||||
|
prod=f2*array_coefs_A(k1,1)*bigI(lambda,mu,l,m,k1,k2,k3)
|
||||||
if (prod == 0.d0) then
|
if (prod == 0.d0) then
|
||||||
cycle
|
cycle
|
||||||
endif
|
endif
|
||||||
@ -2192,6 +2265,7 @@ else
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
ylm_orb=factor*accu
|
ylm_orb=factor*accu
|
||||||
|
deallocate (array_coefs_A)
|
||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user