      program IPT_REAL_AXIS
C======================================================================
C     
C             Hubbard Model on the Bethe Lattice 
C              with infinite Coordination Number
C     
C              - 2 nd order pertubation theory -
C                  away from half filling
C                    MULTI BAND VERSION 
C     (use only for the single-band case because the correlation term
C     appearing in the multiband term is poorly approximated here and
C     should better be approximated using the CPA)
C                  FOR FINITE TEMPERATURE
C                USING dop_mu0t.in FROM T=0 CODE
C
C======================================================================
      
C     ------- set parameters -------
      
      real*8 pi
      parameter(pi=3.141592653589793d0)
      
      integer n,Npoints
      parameter(n=2048*4)
C     (n must be a power of 2.)
      parameter(Npoints=4096)
      
      double complex zi
      parameter(zi=(0.d0, 1.d0))
      
C     ------- model parameters -------
      
      real*8 t,u,umin,umax, ustep
      integer l                 ! Band degeneracy
      integer lh
      
C     ------- Input parameters -------

c      character p*20
      
      integer wnmax, wnmin, stepsi,mi
      integer nloop,inguess,nu
      integer nmu0t             ! Number of mu0 data points read in.
      integer imse_mode
      
C     ------- Dummy variables -------
      
c      character outfile*9
      integer i,i1,i2,k,iloop,index,idx
      real*8 tmp,wout,rout,iout
      double complex ctmp

C     ------- parameters for evaluation -------
      
      integer print_test,stepm,matout,imin

      real*8 Temp,Tmin,Tmax,Tstep
      real*8 dop,dopmin,dopmax,dopstep,dopin
      real*8 mu0t,mu,deltamu,dmuup,dmulow,dmunew,dmuinput
      real*8 deltamu1,deltamu2,deltamu3,deltamu4
      real*8 ynew,yup,ylow,accy,accmu,om(n)
      real*8 average,helpterm2,nb,nb0,corrterm,noccg      
      real*8 nocc,norm,norm2,nocc2
      real*8 domega,tomega, dtime,ttime,time(n),omega(n)
      real*8 acc,eta,nf(n),expon,ommax,matTemp
      real*8 mu0tarray(500),doparray(500),mu0t2array(500)
      real*8 mu0tin,y,theta,helpterm
      real*8 f
      real*8 gss1,gss2

      double complex z(n)
      

C     ------- Green's functions -------
C     ------- index 1 (2) corresponds to '+' ('-') -------

      double complex gtin(2,2,n),gfin(2,2,n)
      double complex gaf1(n),gaf2(n),gaf3(n),gaf4(n)
      double complex gaf(n),gafin(n),gafinvhelp(n)
      double complex g0a(n),g0at(n),g0af(n)
      double complex g0t(2,2,n),gm(n),g0m(n)
      double complex At(n),Af(n),Atf(n)
      double complex Att(n),Attt(n),Attf(n)       
      double complex gfst
      
      
C     ------- Self energy -------

      double complex sigt(2,2,n),sigm(n) 
      double complex sigat(n),sigaf(n),sig0af(n)
      
C     ------- DOS -------

      real*8 dos1(n),dos2(n), dos3(n), dos4(n)
      
      
C     ------- Declaration of global variables -------

      real*8 wfftc(2*n+15), cpy(1*n)
      real*8 aif0(n/2), arf0(n/2), aif1(n/2), arf1(n/2)
      real*8 aif2(n/2), arf2(n/2), aif3(n/2), arf3(n/2)
      real*8 aaf0(n/2), corfaf(n/2)
      real*8 ait0(n/2), art0(n/2), ait1(n/2), art1(n/2)
      real*8 ait2(n/2), art2(n/2), ait3(n/2), art3(n/2)
      real*8 aat0(n/2), corfat(n/2)
      real*8 a01,a0r,a1i,a1r,a2i,a2r,a3i,a3r,aa0,s,c

      common/corre1/ wfftc, cpy 
      common/corre2/ aif0, arf0, aif1, arf1,
     *               aif2, arf2, aif3, arf3,
     *               aaf0, corfaf
      common/corre3/ ait0, art0, ait1, art1,
     *               ait2, art2, ait3, art3,
     *               aat0, corfat
      common/corre4/ a01,a0r,a1i,a1r,a2i,a2r,a3i,a3r,aa0,s,c

C     ------- List of subroutines -------
  
      external init, init_GF
      external fft_ttof, fft_ftot, fftMttof, fftMftot
      external invert
      external writeda
 

C     ====================================================================

c      call getarg(1,p) 
c      decode (20,'(i8)',p) idx
	idx=100

C     ------- Enter Data -------

      open(10, file='input.dat',status='unknown')

      print*, 'Reading data: umin, umax,ustep'
      read(10,*) umin,umax,ustep
      print*, real(umin), real(umax), real(ustep)
      
      print*, 'Reading data: Tmin, Tmax, Tstep'
      read(10,*) Tmin, Tmax, Tstep
      print*, real(Tmin), real(Tmax), real(Tstep)

      print*, 'Reading data: dopmin, dopmax, dopstep'
      read(10,*) dopmin, dopmax, dopstep
      print*, real(dopmin), real(dopmax), real(dopstep)

      print*, 'Reading data: nmu0t'
      read(10,*) nmu0t
	l=1
      print*, l,nmu0t

      print*, 'Reading data: t, nloop'
      read(10,*) t, nloop
      print*, real(t), nloop

      print*, 'Reading data: eta, tomega'
      read(10,*)  eta, tomega
      print*, real(eta), real(tomega)

      print*, 'Enter: accuracy, 0-level,accy,accmu'
      read(10,*) acc,accy,accmu
      print*, real(acc),real(accy),real(accmu)

      print*, 'Enter output data: wnmin, wnmax, stepsi'
      read(10,*) wnmin, wnmax, stepsi
      print*, wnmin, wnmax, stepsi

      print*, 'Output of Matsubara GF:'
      read(10,*) ommax,stepm
      
      print*, 'Enter: initial guess (1=metal, 2=insulator), average'
      read(10,*) inguess, average
      print*, inguess,average

      print*, 
     *  'Print file.test (1:on, 0:off), matout'
      read(10,*) print_test,matout
      print*, print_test,matout

      print*, 'imse_mode (1:on, 0:off)'
      read(10,*) imse_mode
      print*, imse_mode

C     ------- Main program -------

C     ------- initialize parameters etc. -------
      
      call init(domega,tomega,dtime,z,time,omega,ttime,acc,eta)
      print*,'domega=',domega
      print*,'dtime=',dtime
      print*,'domega*dtime=',domega*dtime

      open(7919,file='dop_mu0t.in',status='unknown') 
      do i=1,nmu0t
         read(7919,*) mu0tin,dopin
         mu0tarray(nmu0t+1-i) = mu0tin
         doparray(nmu0t+1-i) = dopin
      enddo
      close(7919)

      write(*,*) "mu0tin=", mu0tin
      write(*,*) "dopin= ", dopin

      call spline(doparray, mu0tarray,nmu0t, 2d30,2d30,mu0t2array)
      do i=1,nmu0t
         mu0t2array(i) = mu0tarray(i)
      enddo

      
C     ------- Test of the spline -------
      
      if (print_test.eq.1) then
         do dop=0.d0,dopin,0.025
            call splint(doparray,mu0tarray,mu0t2array,nmu0t,dop,mu0t)
            write(225,*) real(dop), real(mu0t)
         enddo
      endif
      
C     ----- Here starts the temperature loop -------
      
      do Temp=Tmin, Tmax, Tstep
         
C     ----- Here starts the U loop -------
         
         do u=umin,umax,ustep
            
C     ----- Here starts the doping loop -------
            
            do dop=dopmin,dopmax,dopstep
               
               write(0,*) 'T   =',real(Temp)
               write(0,*) 'U   =',real(u) 
               write(0,*) 'dop =',real(dop)
               
C     ------- determine mu0t, initialize mu -------
               
               call splint(doparray, mu0tarray,mu0t2array,
     *              nmu0t,dop,mu0t)

               mu=mu0t+dfloat(2*l-1)/dfloat(2*l)*u
               
               nocc=(1.d0-dop)/dfloat(2*l)
               deltamu=0.d0

C     ------- Initialize the GF -------
               
               call init_GF(Temp,gtin,gfin,gafin,t,omega,z,
     *                      domega,mu0t,dtime,inguess,u,l)
               
C     ------- Set GF on initial value -------
               do i=1,n
                  gfst=gafin(i)
                  gaf(i)=gfst
                  gaf1(i)=gfst
                  gaf2(i)=gfst
                  gaf3(i)=gfst
                  gaf4(i)=gfst          
               enddo
               deltamu1=0.d0
               deltamu2=0.d0
               deltamu3=0.d0
               deltamu4=0.d0
               
C     ------- Determine the fermi factors -------

               do k=1,n
                  expon=omega(k)/Temp
                  if (dabs(expon).ge.1.d2) then 
                     if (k.eq.n/2) then 
                        nf(k)=0.5d0
                     endif
                     if (k.lt.n/2) then
                        nf(k)=1.d0
                     endif
                     if (k.gt.n/2) then
                        nf(k)=0.d0
                     endif
                  else
                     nf(k)=1.d0/(1.d0+dexp(expon))
                  endif
               enddo
               
C     ------- Here starts the iteration -------
              
               iloop=1
               
 100           print*, (nloop-iloop+1), ' loops remain'
               
C     ------- Normalize G, determine nocc -------
               
               noccg=0.d0
               norm=0.d0
               do k=4,n-2,2
                  tmp = 1.d0/pi*dimag(gaf(k))*domega*2.d0/3.d0
                  noccg = noccg + tmp*nf(k)
                  norm = norm + tmp
               enddo
               do k=3,n-1,2
                  tmp = 1.d0/pi*dimag(gaf(k))*domega*4.d0/3.d0
                  noccg = noccg + tmp*nf(k)
                  norm = norm + tmp 
               enddo
               do k=2,n,(n-2)
                  tmp = 1.d0/pi*dimag(gaf(k))*domega*1.d0/3.d0
                  noccg = noccg + tmp*nf(k)
                  norm = norm + tmp
               enddo
               noccg=noccg/norm

               do k=1,n
                  gaf(k)=gaf(k)/norm
               enddo
               
C     ------- Determine g0a  and then Af,Atf -------

               do k=1,n
                  g0a(k)=1.d0/(z(k)+mu0t-t*t*gaf(k))
                  if (dimag(g0a(k)).lt.0) then
                     g0a(k)=dble(g0a(k))
                  endif
               enddo

               do k=1,n
                  y=1.d0/pi*dimag(g0a(k))
                  Af(k)=y
                  Atf(k)=Af(k)*nf(k)
                  Attf(k)=Af(k)*(1-nf(k))
               enddo
               
               nocc2=0.d0
               norm2=0.d0
               do k=4,n-2,2
                  tmp = 1.d0/pi*dimag(g0a(k))*domega*2.d0/3.d0
                  nocc2 = nocc2 + tmp*nf(k)
                  norm2 = norm2 + tmp
               enddo
               do k=3,n-1,2
                  tmp = 1.d0/pi*dimag(g0a(k))*domega*4.d0/3.d0
                  nocc2 = nocc2 + tmp*nf(k)
                  norm2 = norm2 + tmp
               enddo
               do k=2,n,(n-2)
                  tmp = 1.d0/pi*dimag(g0a(k))*domega*1.d0/3.d0
                  nocc2 = nocc2 + tmp*nf(k)
                  norm2 = norm2 + tmp
               enddo
               
               nocc2=nocc2/norm2 
               print*,'nocc(G)=', nocc,'  norm=', norm
               print*,'nocc(G0)=', nocc2,'  norm2=', norm2   
               print*,'delta=', real(dop)         
               
C     ------- Fouriertransform Af, Atf -------

               call fft_ftot(Af, At,domega)
               call fft_ftot(Atf,Att,domega)
               call fft_ftot(Attf,Attt,domega)
               
C     ------- Calculate g0t (matrix!) -------

               do i=1, n
                  if (i.gt.(n/2)) then 
                     theta=1.d0 
                  else 
                     theta=0.d0 
                  endif
                  if (i.eq.n/2) then
                     theta=0.5
                  endif
                  g0t(2,2,i)=-2.d0*pi*zi*
     *                 (theta*Attt(i)-(1.d0-theta)*Att(i))
                  g0t(2,1,i)= 2.d0*pi*zi*Att(i)
                  g0t(1,2,i)= -2.d0*pi*zi*(At(i)-Att(i))
                  g0t(1,1,i)=g0t(1,2,i)+g0t(2,1,i)-g0t(2,2,i)
               enddo
               
C     ------- Calculate the 2 nd order diagram for the self-energy -------

               do i=1,n-1
                  mi=n-i
                  sigt(2,2,i)=u*u*g0t(2,2,i)*g0t(2,2,i)*g0t(2,2,mi)
                  sigt(1,2,i)=-u*u*g0t(1,2,i)*g0t(1,2,i)*g0t(2,1,mi)
               enddo
               sigt(2,2,n)=sigt(2,2,n-1)
               sigt(1,2,n)=sigt(1,2,n-1)
               
               do i=1,n
                  sigat(i)=sigt(2,2,i)+sigt(1,2,i)
                  g0at(i)= g0t(2,2,i)-g0t(1,2,i)
               enddo
               
               sigat(n/2)=dble(sigat(n/2))+0.5*dimag(sigat(n/2-1))*zi
               
C     ------- Fouriertransform g0t,sigt -------
               
               call fft_ttof(sigat,sig0af,dtime)
               
               do i=1,n
                  if (dimag(sig0af(i)).lt.0) then
                     sig0af(i)=dble(sig0af(i))
                  endif
               enddo
               
C     ------- Determine the interpolating self energy -------
               
               nb=nocc
               nb0=nocc2
               lh=2*l-1
               corrterm=0.d0
               if ((dabs(u).gt.1d-5).and.(l.gt.1)) then
                  do k=1,n
                     if (nloop.eq.1) sigaf(k)=sig0af(k)
                     corrterm = corrterm +
     *                    dimag(gaf(k)*(sigaf(k)+U*nb*(2*l-1)))
                  enddo
                  corrterm = corrterm*domega/U/pi*dfloat(lh-1)
C     -------
                  print*,'corrterm=', corrterm
               endif
               helpterm2=(dfloat(lh)*nb-dfloat(lh*lh)*nb*nb+
     *              corrterm)/dfloat(lh)/nb0/(1.d0-nb0)
               if (dabs(u).gt.1d-5) then
                  helpterm=((1.d0-dfloat(lh)*nb)*U-mu+mu0t)/
     *                       dfloat(lh)/nb0/(1.d0-nb0)/u/u
               else
                  helpterm=0.d0
               endif 

C     -------
               print*, 'helpterm=', helpterm, ' helpterm2=', helpterm2
               do i=1,n
                  sigaf(i)=helpterm2*sig0af(i)/(1.d0-helpterm*sig0af(i))
               enddo
               
C     ------- Determine gafinvhelp -------
               
               do i=1,n
                  gafinvhelp(i)=1.d0/g0a(i)-sigaf(i)
               enddo
               
C     ------- determine new mu -------
               
               dmuup = deltamu+U/2.d0
               dmulow = deltamu-U/2.d0
               dmuinput = dmuup
               yup = f(dmuinput,gafinvhelp,nocc,domega, omega,Temp)
               dmuinput = dmulow
               ylow = f(dmuinput,gafinvhelp,nocc,domega,omega,Temp)
               ynew=1.d0
 200           if ((dabs(ynew).gt.accy).and.
     *              (dabs(dmuup-dmulow).gt.accmu)) then
                  dmunew=dmuup-(dmuup-dmulow)/(yup-ylow)*yup
                  ynew=f(dmunew,gafinvhelp,nocc,domega,omega,Temp)
                  if (ynew.gt.0) then 
                     dmuup=dmunew
                     yup=ynew
                  else
                     dmulow=dmunew
                     ylow=ynew
                  endif        
                  goto 200
               endif
               deltamu=dmunew
               
               
               deltamu=(deltamu4*average+deltamu)/(1.d0+average)
               deltamu1=deltamu2
               deltamu2=deltamu3
               deltamu3=deltamu4
               deltamu4=(deltamu1+deltamu2+deltamu3+deltamu)/4.d0
               deltamu=deltamu4
               
               
               mu=mu0t+deltamu+u*nocc*lh  !The factor lh was missing and
C                                         !introduced by E.L. on 11/19/98
               print*, 'mu0t=', mu0t
               print*, 'mu=', mu, '        deltamu=', deltamu
c               open(20000,file='chem.dat',access='append',
c     &              status='unknown')
c               write(20000,*) iloop, mu
c               close(20000)
               
C     ------- Determine new gaf -------
               
               do i=1,n
                  gaf(i)=1.d0/( 1.d0/g0a(i)-sigaf(i)+deltamu )
               enddo
               
C     ------- Reshuffle old solutions and check convergence -------
              
               do i=1,n
                  do i1=1,2
                     do i2=1,2
                        gaf1(i)=gaf2(i)
                        gaf2(i)=gaf3(i)
                        gaf3(i)=gaf4(i)
                        gaf4(i)=(gaf1(i)+gaf2(i)+ 
     *                           gaf3(i)+gaf(i))/4.d0
                        gaf(i)=gaf4(i)               
                     enddo
                  enddo
               enddo
               
               do k=1,n
                  dos1(k)=dos2(k)
                  dos2(k)=dos3(k)
                  dos3(k)=dos4(k)
                  dos4(k)=1.d0/pi*dimag(gaf(k))
               enddo
               
               gss1=0.d0
               gss2=0.d0
               do i=1,n,n/64
                  gss1=gss1+(dos4(i)-dos1(i))**2
               enddo
               if (inguess.eq.1) then      
                  do i=(n/2-200), (n/2+200),10
                     gss2=gss2+(dos4(i)-dos1(i))**2
                  enddo
               else
                  nu=u/2.d0/domega+n/2
                  do i=(nu-200), (nu+200),10
                     gss2=gss2+(dos4(i)-dos1(i))**2
                  enddo 
               endif     
               
               gss1=gss1/64.d0
               gss2=gss2/40.d0
               print*, 'overall accuracy:',dsqrt(gss1)
               if (inguess.eq.1) then
                  print*, 'accuracy around zero:',dsqrt(gss2)
               else
                  print*, 'accuracy around U/2:', dsqrt(gss2)
               endif
               if (((gss1.lt.acc).and.(gss2.lt.acc)).and.
     *              (iloop.gt.6)) then 
                  print*,'Convergence after', iloop, ' iterations'
               else
                  if (iloop.eq.nloop) then 
                     print*, 'no convergence'
                  else
                     iloop=iloop+1
                     goto 100
                  endif
               endif
               
c               open(83,file='nloop.dat',access='append',
c     *                                  status='unknown')
c               write(83,*) u,Temp,mu0t, '-->', iloop
c               close(83)

C     ------- Determine the Matsubara frequency dependent quantities -------
               
               if (matout.eq.1) then
                  matTemp=Temp
                  if (Temp.lt.0.0009) then 
                     matTemp=0.005
                  endif
                  if (ommax.gt.1d10) then
                     imin=1
                  else
                     imin=n/2
                  endif 
                  do i=imin,n,stepm
                     om(i)=dfloat(2*(i-n/2)+1)*matTemp*pi
                  enddo
                  do i=imin,n,stepm
                     if (om(i).lt.ommax) then
                        gm(i)=0.d0
                        g0m(i)=0.d0
                        sigm(i)=0.d0
                        do k=1,n
                           ctmp = 1.d0/(om(i)*(0,1)-omega(k))
                           ctmp = ctmp/pi*domega
 
                           g0m(i) = g0m(i) + ctmp*dimag(g0a(k))
                           sigm(i) = sigm(i)+ ctmp*dimag(sigaf(k))
                           gm(i) = gm(i) + ctmp*dos4(k)*pi
                        enddo
                     endif
                  enddo
               endif
               
C     ------- Write out results -------
               
               index = idx+10000
c               write(outfile,'(''green.'',i3)') idx
               open(index,file='green',status='new')

C     Write out the temperature,chemical potential, 
C     doping, mesh size and spin degeneracy

               write(index,'(''#'',5e18.10)') temp,mu,dop,domega,2.0
               do k = wnmin,wnmax
                  wout = omega(k)
                  rout = real(gaf(k))
                  iout = real(dimag(gaf(k)))
                  write(index,'(3e18.10)') wout,rout,iout
               enddo
               close(index)
               idx = idx + 1
               
               if(print_test.eq.1) then
                  call writeda(wnmax, wnmin,stepsi, dos4,domega,
     *                         Temp,u,sigaf,sig0af,gafin,gaf,g0af,
     *                         af,nocc,mu,mu0t,print_test,gm,om,
     *                         ommax,stepm,g0a,sigat,l,imse_mode,
     *                         g0m,imin,sigm)
                  
               endif

C     ------- Here ends the doping loop -------
            enddo
C     ------- Here ends the U loop -------
         enddo
C     ------- Here ends the temperature loop -------
      enddo
      end
      
C     ================================================================
















