From 64ea60c7277af739abeed790ab9d23b032535bff Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Tue, 7 Apr 2015 16:59:42 +0200 Subject: [PATCH] Working Local --- scripts/get_basis.sh | 6 +- src/MonoInts/int.f90 | 83 +--------- src/MonoInts/pot_ao_ints.irp.f | 63 +++++-- src/Properties/need.irp.f | 289 --------------------------------- 4 files changed, 54 insertions(+), 387 deletions(-) delete mode 100644 src/Properties/need.irp.f 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 -! -!