diff --git a/scripts/get_basis.sh b/scripts/get_basis.sh
index 51b0a4f0..b5be03ac 100755
--- a/scripts/get_basis.sh
+++ b/scripts/get_basis.sh
@@ -42,9 +42,9 @@ then
echo "ERROR"
exit 1
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}
diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90
index 640688d0..dc07fa54 100644
--- a/src/MonoInts/int.f90
+++ b/src/MonoInts/int.f90
@@ -550,7 +550,6 @@ double precision,allocatable :: array_R_loc(:,:,:)
double precision,allocatable :: array_coefs(:,:,:,:,:,:)
double precision int_prod_bessel_loc,binom,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)
@@ -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
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
@@ -1727,87 +1727,6 @@ end
RETURN
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)
implicit none
integer n,k
diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f
index 6ad38d6a..7fd24174 100644
--- a/src/MonoInts/pot_ao_ints.irp.f
+++ b/src/MonoInts/pot_ao_ints.irp.f
@@ -8,7 +8,7 @@
double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3)
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
! Important for OpenMP
@@ -38,13 +38,55 @@
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 DEFAULT (NONE) &
!$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 v_k, n_k, dz_k, klocmax) &
+ !$OMP num_A, num_B, Z, c, n_pt_in) &
!$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
!$OMP DO SCHEDULE (guided)
do j = 1, ao_num
@@ -74,18 +116,13 @@
C_center(1) = nucl_coord(k,1)
C_center(2) = nucl_coord(k,2)
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)
-
- print*, A_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)
-
+ c = c - Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center)
+ c = c - Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,C_center)
enddo
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
enddo
enddo
enddo
diff --git a/src/Properties/need.irp.f b/src/Properties/need.irp.f
deleted file mode 100644
index 22cb6a48..00000000
--- a/src/Properties/need.irp.f
+++ /dev/null
@@ -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
-
-!
-!
-!
-!
-!
-!
-!
-!
- 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
-
-
-!
-!
-!
-!
-!
-!
-!
-!
-!
-!
-!
-! gcf
-! gser
-!
-!
-!
- 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
-!
-!
-
-
-!
-!
-!
-!
-!
-!
-!
-!
-!
-!
-!
-! gammp
-!
-!
-!
- 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
-!
-
-!
-!
-!
-!
-!
-!
-!
-!
-!
-!
-!
-!
-! gammp
-!
-!
-!
- 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
-
-!
-!
- 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
-
-!
-!
-!
-!
-!
-!
-!
-!
- 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
-!
-!