C (C) Copr. 12/2000 Alexander Moroz C Wave-scattering.com C http://www.wave-scattering.com/ C wavescattering@yahoo.com C C Last update: January 2001 C C This is OLA, a computer code to calculate the free-space C Green's function which satisfies C Bloch periodicity condition along one direction of the C two-dimensional space. This code is based on the article C A. Moroz, Exponentially convergent lattice sums, C Opt. Lett. 26, 1119-1121 (2001). C C This code is made available on the following conditions: C C - You can use and modify it, as long as acknowledgment is given C to the author (that's me). I'd appreciate a reprint or e-copy C of any paper where my code has been used. C C - You can not sell it or use it in any way to get money out of it C (save for the research grants you might get thanks to it) C C - I'm making this code available to the best of my knowledge. C It has been tested as throughly as possible. C It work for Windows and Linux platforms. However, it's C your responsibility to check for its accuracy in the environment C and size/shape ranges you use. Some compilers C (e.q. Compaq Visual Fortran) do not like the following part of the C DL1F1IN2-routine: C C DO 30 IK=0,NK C if (ik.eq.0) go to 5 C do 25 ij=-1,1,2 C 5 qpk=Q+dble(ij*IK)*2.d0*pi/scale C C which makes a jump from outside into the "do 25" loop. (The "if" here C is to insure that the "do 25" part for ik=0 is only C executed a single time.) C The variable"ij" then becomes a "random" variable and no converged C output is produced. C C C By default, x=0.2d0 and y=0.03d0. C If you like to enter the (x,y) values interactively, C go to routine GFF1IN2 and uncomment the following lines therein: C C c write(6,*)'Read in the x-component of R' C c read(5,*) x C x=0.2d0 C c write(6,*)'Read in the y-component of R' C c read(5,*) y C y=0.03d0 C C (and, of course, comment the x=0.2d0 and y=0.03d0 lines). C C Angular-momentum cut off LMAX: C LMAX is set by default to LMAX=50 in main. The same value of LMAX C has been employed to generate results in my Opt. Lett. 26, 1119-1121 (2001). C Since dependent routines are dimensioned to up to LMAX=50, you can C straightforwardly, and with impunity, vary LMAX in main from 1 up to 50. C Above that cutoff: C C 1) you have to raise the parameter value of lbess in C the routine GFF1IN2 so that C LMAX.LE.LBESS C and C C 2) you have to raise the dimensions of factorial common field C fac(0:50) to at least fac(0:LMAX) and initialize the field C up to LMAX (follow "factorial initialization:" in main) C C - Disclaimer. The usual stuff. I claim no responsibility for any C misuse or damage arising from its use. If this program erases your C hard disk, makes your computer cry out for Linux, has you fired, C makes you lose your boy/girlfriend, melts your TV set, spends out C your lifesaving or forces you to get to decaf coffee, it's your C problem. For your tranquility, no such side effects has been C reported ... yet. C C WHILE THE COMPUTER PROGRAMS HAVE BEEN TESTED FOR A VARIETY OF CASES, C IT IS NOT INCONCEIVABLE THAT THEY MAY CONTAIN UNDETECTED ERRORS. ALSO, C INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF C VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON, C THE AUTHORS DISCLAIM ALL LIABILITY FOR C ANY DAMAGES THAT MAY RESULT FROM THE USE OF THE PROGRAMS. C I would highly appreciate informing me of any problems encountered C with this code. C C Improvements, comments, additions: all welcome at C C wavescattering@yahoo.com C C ======================================================= C C REQUIREMENTS & NOT DISTRIBUTED DEPENDENT ROUTINES: C C 1) You should supply your own gnzcylf(zin,lmax,jl,djl,nl,dnl) C routine which generates the Bessel functions JL(0:lbess),NL(0:lbess) C of the first and second order and their derivatives C DJL(0:lbess),DNL(0:lbess) C for a complex argument ZIN. The "gnzcylf" is to be used in GFF1IN2. C C My programs was tested together with the AMOS package of C complex Bessel functions routines which C you can, together with their dependencies, C freely download at Netlib (http://www.netlib.org/amos). C Before using the AMOS package, machine constants C in the real function D1MACH(I) and in the integer function I1MACH(I) C have to be adjusted to your computer. When using Microsoft Windows, C you can activate the same machine constants as listed therein under C the option of Silicon Graphics IRIS. C C 2) For running of the program you will also need some of the C copyrighted routines from Numerical Recipes Software, which can be found C on http://www.nr.com C (I have not received any reply yet from Numerical Recipes C Software on my request for a permission to distribute the routines.) C The routines have to be transformed into double precision. In the case of C the Numerical Recipes routines C qrombt(func,a,b,ss), polint(xa,ya,n,x,y,dy), trapzd(func,a,b,s,n), C follow the instructions on C http://www.wave-scattering.com/dlsum1in2.html C which also provides some supplementary information. C For the remaining routines (erfc.f with its dependencies), C you can either follow instructions written down in C Numerical Recipes or you can copy and paste the C following declarations: C C DOUBLEPRECISION FUNCTION erfc(x) C REAL*8 x,sqrpi C CU USES gammp,gammq C REAL*8 gammp,gammq C * C DOUBLEPRECISION FUNCTION gammp(a,x) C REAL*8 a,x C CU USES gcf,gser C REAL*8 gammcf,gamser,gln C * C SUBROUTINE gcf(gammcf,a,x,gln) C implicit none C INTEGER ITMAX C REAL*8 a,gammcf,gln,x,EPS,FPMIN C PARAMETER (ITMAX=100,EPS=1.d-15,FPMIN=1.d-30) C CU USES gammln C INTEGER i C REAL*8 an,b,c,d,del,h,gammln C * C SUBROUTINE gser(gamser,a,x,gln) C implicit none C INTEGER ITMAX C REAL*8 a,gamser,gln,x,EPS C PARAMETER (ITMAX=100,EPS=1.d-15) C CU USES gammln C INTEGER n C REAL*8 ap,del,sum,gammln C * C DOUBLEPRECISION FUNCTION gammln(xx) C implicit none C REAL*8 xx C INTEGER j C DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) C * C DOUBLEPRECISION FUNCTION gammq(a,x) C REAL*8 a,x C CU USES gcf,gser C REAL*8 gammcf,gamser,gln C C C ======================================================= PROGRAM OLA C--------/---------/---------/---------/---------/---------/---------/-- C 1D LATTICE ORIENTED ALONG THE X-AXIS C--------/---------/---------/---------/---------/---------/---------/-- IMPLICIT NONE integer lmax,nb,ncmb,nr,nk integer NFIN,NOUT,NSTACK,NTEST real*8 Q0,SCALE,phin logical ynbrug * * number of the output unit PARAMETER (NOUT=60) * total amount of test repetitions PARAMETER (NTEST=1) * angular-momentum cut off PARAMETER(LMAX=50) * the number of atoms within the unit cell PARAMETER(NB=1) * the number of atom pairs within the unit cell PARAMETER(NCMB=1+NB*(NB-1)/2) * cut off on the max. length of the direct lattice vectors PARAMETER(NR=50) * cut off on the max. length of the reciprocal lattice vectors PARAMETER(NK=50) * the length of the direct lattice unit cell PARAMETER(SCALE=1.d0) * the Ewald parameter PARAMETER(Q0=.011d0) * integer il,ilp,ib real*8 q,sg,qscale,pi,fac(0:50) COMPLEX*16 DL1(0:LMAX,NCMB),DL2(0:LMAX,NCMB),B(0:LMAX,NCMB) * common/fac50/fac * DATA PI/3.141592653589793d0/ * do 15 itest=1,ntest * * factorial initialization: * FAC(0)=1.D0 DO 1 il=1,50 FAC(il)=DBLE(il)*FAC(il-1) 1 CONTINUE * C======================================================================= * Setting up the components: * qscale=2.d0*pi/scale ! Basis vector of the dual/reciprocal lattice * * Setting scattering angle and perpendicular momentum component: * sg=qscale/0.23d0 phin=pi/8.d0 q=sg*sin(phin) ! perpendicular momentum component c q=sg*cos(phin) * * selecting the Bloch vector within the 1st Brillouin zone: do il=0,nk if (abs(q).gt.qscale) q=q-qscale if (abs(q).lt.qscale) go to 5 enddo c write(6,*)'Incorrect Bloch momentum' c stop 11 c end if * 5 CONTINUE C======================================================================= * Checking set up: * if(phin.gt.pi/2.d0) then write(6,*)'PHIN=',phin,'.gt.pi/2.d0 in OLA' stop end if * if(sg.gt.nk*.9d0) then write(6,*)'NK=',nk,'.lt.sg in OLA' write(6,*)'Raise NK in OLA so that NK>1.3*sg' stop end if C======================================================================= * Execution * * call BLF1IN2(Q0,Q,SG,SCALE,LMAX,NCMB,NK,NR,DL1,DL2,B) * * call gff1in2(SG,LMAX,NCMB,B) * 15 CONTINUE *--------/---------/---------/---------/---------/---------/---------/-- end *------------------------------------------------------------------------------------- SUBROUTINE BLF1IN2(Q0,Q,SG,SCALE,LMAX,NCMB,NK,NR,DL1,DL2,B) C--------/---------/---------/---------/---------/---------/---------/-- C >>> LMAX,NCMB,DL1,DL2,DL3 C <<< BL C =================================== C g_{l_1-l_2}=i^{|l_1|-|l_2|} \sum_l i^{-|l|} \dt_{l,l_1-l_2} D_l C g_{-l,-l'}=(-1)^{l+l'} g_{ll'}^* C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer lmax,ncmb,nk,nr,il,ib real*8 sg,scale,xtol,xxtol,q,q0,xetest real*8 DL3 COMPLEX*16 DL1(0:LMAX,NCMB),DL2(0:LMAX,NCMB) COMPLEX*16 B(0:LMAX,NCMB) * call dl1f1in2(Q0,Q,SG,SCALE,NK,NCMB,LMAX,DL1) call dl2f1in2(Q0,Q,SG,SCALE,NR,NCMB,LMAX,DL2) call dl3in2(Q0,SG,DL3) * DO 2 IB=1,NCMB DO 1 IL=0,LMAX C--------/---------/---------/---------/---------/---------/---------/-- * security trap for numerical stability of the Ewald summation: * if (dble(DL1(IL,IB))*dble(DL2(IL,IB)).lt.0.) then xetest=1.d0-abs(DL1(IL,IB)/DL2(IL,IB)) if(abs(xetest).lt.4.d-1) then write(6,*)'In BLF1IN2:' write(6,*)'The Ewald summation may be instable - change Q0' end if end if * B(IL,IB)=DL1(IL,IB)+DL2(IL,IB) 1 CONTINUE 2 CONTINUE B(0,1)=B(0,1)+dl3 *--------/---------/---------/---------/---------/---------/---------/-- return end *------------------------------------------------------------------------------------- SUBROUTINE GFF1IN2(SG,LMAX,NCMB,B) C--------/---------/---------/---------/---------/---------/---------/-- C >>> SG,LMAX,NCMB,B C <<< ZGRF C The atan() function returns the arc tangent in radians and C the value is mathematically defined to be between -PI/2 C and PI/2 (inclusive). C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer lbess *::: cutoff on bessel functions PARAMETER (lbess=50) integer lmax,ncmb,il,ib real*8 sg,scale,xtol,xxtol,pi,phi,rvs,x,y cc REAL*8 JL(0:lbess),NL(0:lbess),DJL(0:lbess),DNL(0:lbess) COMPLEX*16 JL(0:lbess),NL(0:lbess),DJL(0:lbess),DNL(0:lbess) COMPLEX*16 B(0:LMAX,NCMB),CI,zyl,zgrf,zin * DATA PI/3.141592653589793d0/ DATA CI/(0.0D0,1.0D0)/ * * checking set up if (lmax.gt.lbess) then write(6,*)'Raise the value of LBESS in GFF1IN2' write(6,*)'Should be LBESS.GE.LMAX' stop end if * c write(6,*)'Read in the x-component of R' c read(5,*) x x=0.2d0 c write(6,*)'Read in the y-component of R' c read(5,*) y y=0.03d0 * rvs=sqrt(x**2+y**2) * cc call gncylf(sg*rvs,lmax,jl,djl,nl,dnl) zin=sg*rvs call gnzcylf(zin,lmax,jl,djl,nl,dnl) * if (x.eq.0.d0) then if(y.gt.0.d0) phi=pi/2.d0 if(y.lt.0.d0) phi=-pi/2.d0 else if ((x.eq.0.d0).and.(y.eq.0.d0)) then write(6,*)'Incorrect input in GFF1in2!!!' stop else if ((x.ne.0.d0).and.(y.eq.0.d0)) then if(x.gt.0.d0) phi=0.d0 if(x.lt.0.d0) phi=pi else if ((x.gt.0.d0).and.(y.ne.0.d0)) then phi=atan(y/x) else if ((x.lt.0.d0).and.(y.ne.0.d0)) then if (y.lt.0.d0) phi=-pi+atan(y/x) if (y.gt.0.d0) phi=pi+atan(y/x) end if * write(6,*)'phi=',phi zgrf=(0.d0,0.d0) * ib=1 * DO 1 IL=0,LMAX zyl=exp(ci*dble(il)*phi) zgrf=zgrf+B(IL,IB)*zyl*jl(il) if (il.ne.0) zgrf=zgrf+B(IL,IB)*dconjg(zyl)*jl(il) 1 CONTINUE zgrf=nl(0)/4.d0+zgrf/sqrt(2.d0*pi) write(6,*)'Green function=',zgrf *--------/---------/---------/---------/---------/---------/---------/-- return end *------------------------------------------------------------------------------------- SUBROUTINE DL1F1IN2(Q0,Q,SG,SCALE,NK,NCMB,LMAX,DL1) C--------/---------/---------/---------/---------/---------/---------/-- C >>> Q0,Q,SG,NK,NB,SCALE,LMAX C <<< DL1 C TOL=1.d-14 C =================================== C 1D LATTICE ORIENTED ALONG THE X-AXIS C Q0 ... Ewald parameter C Q ... (Parallel component of) the Bloch vector C SCALE ... length (volume) of the unit lattice cell C NK ... cutoff on the momentum vectors C NCMB ... the number of atoms pairs within the unit cell C LMAX ... cutoff on the angular momenta C C Q0 determines NKMAX for which zint is not negligible *--------/---------/---------/---------/---------/---------/---------/-- implicit none real*8 TOL character*1 yntest C ::: relative error. If the convergence * within TOL is not reached, program issues warning PARAMETER (TOL=1.d-14) * Should convergence test be performed??? PARAMETER (yntest='n') * integer nk,ncmb,lmax,ij,il,ila,lm2n,ib,ik,in,inp1 real*8 pi,sg,xlf,scale,xtol(0:lmax),xxtol(0:lmax),xint,xpref REAL*8 Q,Q0,qpk,rqpk,erfc,funka,fac(0:50) complex*16 kperp,zint0,zint,zar,zarsq,zpref,zdl1 COMPLEX*16 DL1(0:LMAX,NCMB),CI * common/fac50/fac * external erfc,funka * DATA PI/3.141592653589793d0/ DATA CI/(0.0D0,1.0D0)/ * * Initialization: if ((yntest.eq.'y').or.(yntest.eq.'Y')) then do il=0,lmax xtol(il)=10.d0*tol enddo end if * do 50 ib=1,ncmb * c DL1 initialization: * DO 3 IL=0,LMAX DL1(IL,IB)=dcmplx(0.d0,0.d0) 3 CONTINUE * ij=1 C ***************************************************************** C >>> STARTS reciprocal lattice summation loop <<< C DO 30 IK=0,NK if (ik.eq.0) go to 5 do 25 ij=-1,1,2 5 qpk=Q+dble(ij*IK)*2.d0*pi/scale rqpk=abs(qpk) * if(sg.gt.rqpk) then kperp=sqrt(sg**2-qpk**2) call qrombt(funka,0.d0,abs(kperp)*sqrt(q0/2.d0),xint) zint0=sqrt(pi)+2.d0*ci*xint else if(sg.lt.rqpk) then kperp=ci*sqrt(qpk**2-sg**2) xint=erfc(abs(kperp)*sqrt(q0/2.d0)) zint0=sqrt(pi)*xint else if(sg.eq.rqpk) then kperp=0.d0 zint0=sqrt(pi) end if zar=-q0*kperp**2/2.d0 zarsq=-ci*kperp*sqrt(q0/2.d0) * C ----------------------------------------------------------------- C >>> STARTS L-loop. Assigns values of DL1 for different l C ----------------------------------------------------------------- do 20 il=0,lmax xpref=1.d0 ila=abs(il) zint=zint0 * do 10 in=0,ila/2 lm2n=ila-2*in * spherical harmonics part: !lattice along the x-direction * if (qpk.gt.0) xlf=1.d0 if (qpk.lt.0) xlf=(-1.d0)**lm2n * zdl1=xpref*xlf*zint*kperp**(2*in-1)*rqpk**lm2n dl1(il,ib)=dl1(il,ib)+zdl1 * c write(6,*) 'abs(dl1(il,ib))=',abs(dl1(il,ib)) * if convergence test is on: if ((yntest.eq.'y').or.(yntest.eq.'Y')) then if (abs(dl1(il,ib)).ne.0.) then xxtol(il)=abs(zdl1)/abs(dl1(il,ib)) if (xxtol(il).lt.xtol(il)) xtol(il)=xxtol(il) end if end if * if(in.ge.ila/2) goto 20 * prefactor: inp1=in+1 xpref=fac(ila)/(fac(lm2n-2)*fac(inp1)*4.d0**inp1) * recursion for zint: if (kperp.ne.0.d0) & zint=(zint-zarsq*exp(-zar)/zar**inp1)/dble(.5d0-inp1) *--------/---------/---------/---------/---------/---------/---------/-- 10 continue * 20 continue 25 continue 30 continue ! end of lattice summation loop * if (yntest.eq.'y'.or.yntest.eq.'Y') then do il=0,lmax if (xtol(il).gt.tol) then write(6,*)'Warning! For il=',il,' and ib=',ib write(6,*)'Convergence in DL1f1IN2.F not reached!' write(6,*)'xtol(il)=',xtol(il) stop end if enddo end if * * Adding l-dependent prefactors: zpref=-sg/(sqrt(2.d0)*scale) do il=0,lmax ila=abs(il) zpref=ci*zpref/sg dl1(il,ib)=zpref*dl1(il,ib) enddo * 50 continue * *--------/---------/---------/---------/---------/---------/---------/-- return end * DOUBLEPRECISION FUNCTION funka(x) real*8 x funka=exp(x**2) return END *------------------------------------------------------------------------------------- SUBROUTINE DL2F1IN2(Q0,Q,SG,SCALE,NR,NCMB,LMAX,DL2) C--------/---------/---------/---------/---------/---------/---------/-- C >>> Q0,Q,SG,NR,SCALE,LMAX C <<< DL2 C TOL=1.d-14 C =================================== C 1D LATTICE ORIENTED ALONG THE X-AXIS C Q0 ... Ewald parameter C Q ... (Parallel component of) the Bloch vector C SCALE ... length (volume) of the unit lattice cell C NR ... cutoff on the latice vectors C LMAX ... cutoff on the angular momenta C NCMB ... the number of atoms pairs within the unit cell C For IB >= 2: R(3,NR,IB) = R(3,NR,1) - TAU(IB) C DL2 C * * * * * * * * * C >>> NUMERICAL CONSTANTS: PI C ============================== C INTERNAL FUNCTIONS: C C generic cmplx(*) ---> DCMPLX(*) C >>> specific CDEXP ---> ZEXP C DATA for TOL, ZERO ! C ============================== C--------/---------/---------/---------/---------/---------/---------/-- implicit none real*8 TOL character*1 yntest C ::: relative error. If the convergence * within TOL is not reached, program issues warning PARAMETER (TOL=1.d-14) * Should convergence test be performed??? PARAMETER (yntest='n') * integer nr,ncmb,ij,il,ib,ir,lmax real*8 Q,Q0,qdr,pi,prefl,xx,rs,sg,al,xlf,scale,xint REAL*8 ZLR(0:LMAX),funka0,funka1,xtol(0:lmax),xxtol(0:lmax) COMPLEX*16 DL2(0:LMAX,NCMB),CI,cef * common/DDL2/ xx * external funka0,funka1 * DATA PI/3.141592653589793d0/ DATA CI/(0.0D0,1.0D0)/ * * Initialization: if ((yntest.eq.'y').or.(yntest.eq.'Y')) then do il=0,lmax xtol(il)=10.d0*tol enddo end if * do 40 ib=1,ncmb * c DL2 initialization: do il=0,lmax dl2(il,ib)=dcmplx(0.d0,0.d0) enddo C ***************************************************************** C >>> STARTS direct lattice summation loop <<< C DO 20 IR=1,NR do 20 ij=-1,1,2 QDR=Q*dble(ij*IR)*scale CEF=EXP(-CI*QDR) rs=ir*scale xx=sg*rs al=q0*sg**2/2.d0 * Calculating the integral: * The recurrence initialization terms: * call qrombt(funka0,0.d0,al,xint) zlr(0)=xint call qrombt(funka1,0.d0,al,xint) zlr(1)=xint * * Using the recurrence relation: * do il=1,lmax-1 zlr(il+1)=dble(il)*zlr(il) & -zlr(il-1)+exp(al-xx**2/(4.d0*al))/al**il zlr(il+1)=zlr(il+1)/(xx/2.d0)**2 enddo C ----------------------------------------------------------------- C >>> STARTS L-loop. Assigns values of DL2 for different l C ----------------------------------------------------------------- do 10 il=0,lmax * spherical harmonics part: !lattice along the x-direction * if (ij.gt.0) xlf=1.d0 if (ij.lt.0) xlf=(-1.d0)**il * dl2(il,ib)=dl2(il,ib) + zlr(il)*CEF*xlf*xx**il * if convergence test is on: if ((yntest.eq.'y').or.(yntest.eq.'Y')) then if(abs(dl2(il,ib)).ne.0.) then xxtol(il)=abs(zlr(il)*CEF*xlf*xx**il)/abs(dl2(il,ib)) if (xxtol(il).lt.xtol(il)) xtol(il)=xxtol(il) end if end if * 10 continue 20 continue * if (yntest.eq.'y'.or.yntest.eq.'Y') then do il=0,lmax if (xtol(il).gt.tol) then write(6,*)'Warning! For il=',il,' and ib=',ib write(6,*)'Convergence in DL2f1IN2.F not reached!' write(6,*)'xtol(il)=',xtol(il) stop end if enddo end if * * Adding l-dependent prefactors: prefl=-1.d0/(2.d0*sqrt(2.d0*pi)) do il=0,lmax dl2(il,ib)=prefl*dl2(il,ib) prefl=-1.d0*prefl/2.d0 enddo * 40 continue * *--------/---------/---------/---------/---------/---------/---------/-- return end * DOUBLEPRECISION FUNCTION funka0(x) real*8 x,xx common/DDL2/ xx if (x.lt.1.d-50) then funka0=0.d0 else funka0=exp(x - xx**2/(4.d0*x))/x end if return END * * DOUBLEPRECISION FUNCTION funka1(x) real*8 x,xx common/DDL2/ xx if (x.lt.1.d-50) then funka1=0.d0 else funka1=exp(x - xx**2/(4.d0*x))/x**2 end if return END *------------------------------------------------------------------------------------- subroutine dl3in2(Q0,SG,DL3) C--------/---------/---------/---------/---------/---------/---------/-- C >>> SG,Q0 C <<< DL3,DL3P C =================================== C DL3 EVALUATES THE THIRD CONTRIBUTION TO THE STRUCTURE CONSTANTS C FOR POSITIVE ENERGIES FOR A 1D LATTICE IN TWO DIMENSION C C SG ... THE DIMENSIONLESS ENERGY PARAMETER (E=EPS*(2PI/A)**2) C Q0 ... THE EWALD PARAMETER C DL3 ... THE THIRD CONTRIBUTION TO THE STRUCTURE CONSTANTS C DL3P... ENERGY DERIVATIVE OF THE THIRD CONTRIBUTION TO THE STRUCTURE C CONSTANTS C * * * * * * * * * C >>> NUMERICAL CONSTANTS: PI,GM, TOL=1.D-5 C ============================== C INTERNAL FUNCTIONS: C C GENERIC SQRT,LOG C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer n real*8 TOL C ::: relative error. If the convergence * within TOL is not reached, program issues warning PARAMETER (TOL=1.d-14) real*8 pi,al,gm,xx,dl3,sg,q0 *--------------------------------------------------------------- c Data: DATA PI/3.141592653589793d0/ DATA GM/5.772156649015328d-1/ !The Euler gamma constant *oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo * Execution al=q0*sg**2/2.d0 xx=al dl3=xx do n=1,500 xx=dble(n)*xx*al/dble(n+1)**2 dl3=dl3+xx if (xx/dl3.lt.tol) go to 5 enddo write(6,*)'Warning! Convergence in DL3IN2.F not reached!' stop 5 continue dl3=gm+log(al)+dl3 dl3=-dl3/(2.d0*sqrt(2.d0*pi)) *--------/---------/---------/---------/---------/---------/---------/-- return end C (C) Copr. 12/2000 Alexander Moroz