c subprograms to calculate HSM-model s-wave scattering for B phase c at all temperatures. c interactive main program in file mhsm.f c This is an important subprogarm in the loci and nhsm program c KUTSU ALIOHJELMAA bgap JOSTAIN OMASTA OHJELMASTASI. c ARGUMENTTEINA OVAT: c sisaan menevat c LAMPOTILA temp=T/Tc, c EPAPUHTAUSPERAMETRIT scatt (ksi0/ltr) ja sigmat (born/unit.) c tarkkuusparametri accur c nasymp=0: no asymptotic corrections to the matsubara sums c =1: asymptotic 1/eps**2 form of matsubara sums taken into account c =else: asymptotic 1/eps**2 and 1/eps**3 forms taken into account c see (ajop5.25 and 30) or (13.116) c ulos tulevat c GAP delta c VEKTORI d, JOSSA TERMIT ind(em) c SUMMAUSRAJA nfreq c transition temperature tc c all real quantites in units of T_c0 subroutine gapimp(temp,scatt,sigmat,accur,nasymp,gap, x d,nfreq,tc) implicit none integer iq, nfreq,idime,iter,nasymp,nsplit,nsplir parameter(idime=52) double precision dold, delta, r, s, temp, em, x(idime) double precision pi, gap,small, gapf, dgapf, ff, ffd double precision sum, sigmat, scatt, d(idime),dt(idime) double precision accur,tc,piit2,gapt parameter(pi = 3.141592654) parameter(small = 0.00001) piit2=2*pi*temp call cutoff(temp,accur,nfreq) if(nfreq.gt.idime)then write(6,*)' too large nfreq *******************************' nfreq=idime end if write(6,3)nfreq,accur,temp 3 format(' nfreq=',i3,' for accur=',f12.8,' temp=',f8.6) call tcimp(scatt,tc) write(6,5)tc,scatt 5 format(' tc=',f8.6,' for scatt=',f8.6) nsplir=max(nfreq,10) if(temp.ge.tc)then gap=0. do 30 iq=1,nsplir d(iq)=pi*scatt 30 continue write(6,2)gap,d(1),d(2),d(3),nsplir return end if c Risto's subprogram call gap_b(nsplir,temp,gap,x,scatt,sigmat) do 40 iq = 1,nsplir d(iq) = x(iq) - piit2*(float(iq)-0.5) 40 continue write(6,2)gap,d(1),d(2),d(3),nsplir 2 format('Risto:gap=',f8.6,' d(1..3)=',3f8.6,' nsplit=',i3) c Tero's subprogram nsplit=nfreq iter=0 r = 1. - sigmat s=scatt/(2.*temp) c normaalitilan alkuarvot do 50 iq=1,nsplit x(iq)=float(iq)-0.5+s dt(iq)=piit2*s 50 continue delta = 1.7638*sqrt(tc - temp)/piit2 c call init (nsplit, delta,r,s,temp,x,idime) 60 iter=iter+1 call newton (nsplit,temp, delta, r, s, small, sum, x,idime) c NEWTONIN ITERAATIO GAPILLE if(nasymp.eq.0)then gapf =log(temp)+sum else if(nasymp.eq.1)then gapf =log(temp)+sum+s/float(nsplit) else gapf =log(temp)+sum+s/float(nsplit) x +(delta**2-2.*s**2)/float(4*nsplit*2) end if dgapf = 0. do 100 iq = 1,nsplit dgapf=dgapf+delta*(1/sqrt(x(iq)**2+delta**2))**3 100 continue dold = delta delta = dold - gapf/dgapf c write(6,4)delta,sum,x(0),x(1),x(2) c4 format('gap,sum,xi(0..2)',5f10.5) if (abs(delta - dold) .gt. small) then goto 60 else c LASKETAAN TERMIT ind(em) do 200 iq = 1,nsplit dt(iq) = piit2*(x(iq) - float(iq)+0.5) 200 continue gapt=piit2*delta write(6,1)gapt,dt(1),dt(2),dt(3),nsplit,iter 1 format('Tero: gap=',f8.6,' d(1..3)=',3f8.6,' nsplit=',i3, x ' iter=',i3) return end if end double precision function ff (x, em, delta, r, s) implicit none double precision em, delta, r, s, x c x:N MAARAAVA YHTALO ff = (x - em)*(r*delta*delta + x*x) - s*x* % sqrt(x*x + delta*delta) end double precision function ffd (x, em, delta, r, s) implicit none double precision em, delta, r, s, x c x:N MAARAAVAN YHTALON DERIVAATTA ffd = r*delta*delta + x*x + (x - em)*2.*x - s* % (2.*x*x + delta*delta)/sqrt(x*x + delta*delta) end subroutine newton (nfreq, t, delta, r, s, small, sum, x,idime) implicit none integer nfreq, iq,idime double precision t, delta, r, s, small, x(idime) double precision em, prev, psum, sum, ff, ffd, pi sum = 0. do 270 iq = 1,nfreq em = float(iq)-0.5 260 prev = x(iq) c TARKENNETAAN HAARUKOITUJA x(iq):N ARVOJA ITEROIMALLA x(iq)=prev-ff(prev,em,delta,r,s)/ffd(prev,em,delta,r,s) if (abs(x(iq) - prev).gt.0.01*small) then goto 260 end if psum = 1./em-1./sqrt(x(iq)**2 + delta**2) sum = sum + psum 270 continue return end subroutine init (nfreq, delta, r, s, t, x,idime) c nykyisellŠŠn tŠmŠ aliohjelma ei ole kŠytšssŠ implicit none integer N, i, iq, nfreq,idime double precision TOL, delta, em, r, s, x(idime) double precision t, p, a, b, gap, fprime, err double precision ff,pi c N ON MAKSIMI HAARUKOINTI M€€R€ JA TOL HAARUKOINTI TARKKUUS N = 30 TOL = 0.01 c ALUSTETAAN TAULUKKO x HAARUKOIMALLA do 400 iq = 1,nfreq a = dble(-100000) b = dble(100000) em = float(iq)-0.5 i = 1 300 p = a + (b - a)/2.0 if ((ff(p,em,delta,r,s).eq.0.).or.(abs((b-a)/2.).lt.TOL)) then x(iq) = p goto 400 else i = i + 1 if (ff(p,em,delta,r,s)*ff(a,em,delta,r,s) .gt. 0) then a = p else b = p end if end if if (i .le. N) then goto 300 end if 400 continue return end subroutine suprg(temp,scatt0,sigmat,bulkg,bulkr,nfreq, x idime,decay1,decay2) c the integrated suppression of the gap around one impurity (ajop5.29) c if the scattering were smeared out, mostly obsolete integer idime,nfreq,ie parameter(pii=3.1415926535897932) double precision temp,scatt0,sigmat,bulkg,bulkr(idime) double precision supint,e0,sq,piit2,sum,xs,decay1,decay2 piit2=2.*pii*temp fprime=0. do 30 ie=1,nfreq e0=float(ie)-0.5 sq=sqrt((piit2*e0+bulkr(ie))**2+bulkg**2) fprime=fprime+piit2*bulkg/sq**3 30 continue sum=0. do 60 ie=1,nfreq e0=float(ie)-0.5 xs=(piit2*e0+bulkr(ie))**2 sum=sum+xs/((xs+(1.-sigmat)*bulkg**2)*(xs+bulkg**2)) 60 continue decay1=pii/sqrt((piit2*(float(nfreq)-0.5))**2+bulkg**2) decay2=pii/sqrt((piit2*(float(nfreq)+0.5))**2+bulkg**2) return end