diff --git a/src/nuclei/nuclei.irp.f b/src/nuclei/nuclei.irp.f index 39e6b399..523bd80f 100644 --- a/src/nuclei/nuclei.irp.f +++ b/src/nuclei/nuclei.irp.f @@ -42,6 +42,7 @@ BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num,3) ] ' Atom Charge X Y Z ' write(6,ft) & '================','============','============','============','============' + do i=1,nucl_num write(6,f) nucl_label(i), nucl_charge(i), & nucl_coord(i,1)*a0, & @@ -52,6 +53,19 @@ BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num,3) ] '================','============','============','============','============' write(6,'(A)') '' + double precision :: dist_min, x, y, z + dist_min = huge(1.d0) + do i=1,nucl_num + do j=i+1,nucl_num + x = nucl_coord(i,1)-nucl_coord(j,1) + y = nucl_coord(i,2)-nucl_coord(j,2) + z = nucl_coord(i,3)-nucl_coord(j,3) + dist_min = min(x*x + y*y + z*z, dist_min) + enddo + enddo + write(6,'(A,F12.4,A)') 'Minimal interatomic distance found: ', & + dsqrt(dist_min)*a0,' Angstrom' + endif IRP_IF MPI_DEBUG diff --git a/src/utils/need.irp.f b/src/utils/need.irp.f deleted file mode 100644 index 0e36351c..00000000 --- a/src/utils/need.irp.f +++ /dev/null @@ -1,210 +0,0 @@ - - double precision function SABpartial(zA,zB,A,B,nA,nB,gamA,gamB,l) - implicit double precision(a-h,o-z) - dimension nA(3),nB(3) - dimension A(3),B(3) - gamtot=gamA+gamB - SABpartial=1.d0 - - 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.include.F' - 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 none - double precision, intent(in) :: x - double precision, external :: gammp - if(x < 0.d0)then - erf0=-gammp(0.5d0,x*x) - else - erf0=gammp(0.5d0,x*x) - endif - end - - - double precision function gammp(a,x) - implicit none - double precision, intent(in) :: a, x - double precision :: gln, gammcf - if(x<0.d0.or.a<=0.d0) then - stop 'error in gammp' - endif - if(x < a+1.d0)then - call gser(gammp,a,x,gln) - else - call gcf(gammcf,a,x,gln) - gammp=1.d0-gammcf - endif - 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) stop 'error in gser' - 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(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 - - 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=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.d0)then - fac=1.d0/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