program heiswolff5
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     Monte Carlo simulation of clean 3d Heisenberg
! 
!     periodic boundary conditions, sequential update of sites
!     uses LFSR113 random number generator
!     -------------------------------------------------------------------
!     History:
!     
!     heiswolff :  8  May 2003   first version, based on 5d Ising program clean_1
!     heiswolff2:  12 May 2003   introduce site dilution
!     heiswolff3:  13 May 2003   changed order of space and time to be identical to Rastko, 
!                                 fixed several bugs 
!     heiswolff4:  23 May 2003   implemented [log<m>] and correlation length, added presweeps in equilibration
!     heiswolff5:  28 May 2003   switch between canonical and grandcanonical impurities
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Parameters
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      implicit none
      integer,parameter         :: r8b= kind(1.D200)
      integer,parameter         :: i4b= kind(2147483647)
      integer,parameter         :: ilog=kind(.true.)
      real(r8b), parameter      :: pi=3.141592653589D0 

      integer(i4b),parameter    :: L=18, Lt=25, L3=Lt*L*L        ! L=linear system size
      real(r8b), parameter      :: qspace=2*pi/L,qtime=2*pi/Lt   ! minimum q values for correlation length
      integer(i4b),parameter    :: NEQ=100,NMESS=100     ! Monte Carlo sweeps
      real(r8b),   parameter    :: TMAX=1.2075D0, TMIN=1.184999D0, DT=-0.0025D0      ! temperature control
      integer(i4b),parameter    :: NTEMP= 1+(TMIN-TMAX)/DT

      real(r8b),   parameter    :: IMPCONC=0.197
      integer(i4b),parameter    :: NCONF=1000                             ! number of disorder configs 
      logical(ilog), parameter  :: CANON_DIS=.true. 

      integer(i4b),parameter    :: IRINIT=1, IM=2147483647                 ! LFSR stuff

      integer(i4b),parameter    :: NCLUSTER0=60                          ! initial cluster flips per sweep    

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Variable declarations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      real (r8b)      :: sx(0:L3-1),sy(0:L3-1),sz(0:L3-1)  ! Heisenberg spin
      logical(ilog)   :: occu(0:L3-1)                      ! occupation of site with spin 
      real (r8b)      :: sxnew, synew, sznew, slen
      real (r8b)      :: nsx,nsy,nsz, dE                   ! sum over neighboring spins

      real(r8b)       :: mx,my,mz,sweepmag,sweepen                  ! magnetization vector 
      real(r8b)       :: sweepmagqt,sweepmagq2, sweepmagq3,magqt,magqs     ! mag(qtime),mag(qspace), 
      real(r8b)       :: Gt,Gs,Gtcon,Gscon                         ! G(qspace), G(qtime), connected versions
      real(r8b)       :: mag, mag2, mag4, bin, susc,en, en2,sph     ! magnetization, its square, energy
      real(r8b)       :: xit,xis,xitcon,xiscon                      ! correlation lengths in space an time, connected versions
      real(r8b)       :: glxit,glxis,glxitcon,glxiscon              ! global correlation lengths in space an time, connected versions

      real(r8b)       :: confmag(NTEMP),confmag2(NTEMP)
      real(r8b)       :: conf2mag(NTEMP), confmag4(NTEMP)              ! configuration averages
      real(r8b)       :: conflogmag(NTEMP)
      real(r8b)       :: confsusc(NTEMP), confbin(NTEMP)
      real(r8b)       :: conf2susc(NTEMP), conf2bin(NTEMP)              ! conf. av. of squared observables
      real(r8b)       :: confen(NTEMP),confsph(NTEMP)
      real(r8b)       :: conf2en(NTEMP),conf2sph(NTEMP)
      real(r8b)       :: confGt(NTEMP),confGs(NTEMP),confGtcon(NTEMP),confGscon(NTEMP)
      real(r8b)       :: confxit(NTEMP),confxis(NTEMP) 
      real(r8b)       :: confxitcon(NTEMP),confxiscon(NTEMP) 

      integer(i4b)    :: m1(0:L3-1)            ! neighbor table
      integer(i4b)    :: m2(0:L3-1)
      integer(i4b)    :: m3(0:L3-1)
      integer(i4b)    :: m4(0:L3-1)
      integer(i4b)    :: m5(0:L3-1)
      integer(i4b)    :: m6(0:L3-1)

      integer(i4b)    :: iconf,init            ! current disorder config
      integer(i4b)    :: N_IMPSITE,iimp        ! number of impurity sites 

      integer(i4b)    :: i1,i2,i3             ! coordinates
      integer(i4b)    :: is                 ! site index
      integer(i4b)    :: ispace,itime 

      integer(i4b)    :: isweep                ! Monte Carlo sweep
      real(r8b)       :: T, beta               ! Monte Carlo temperature
      integer(i4b)    :: itemp

      integer(i4b)    :: nsweepflip,ncluster
      real(r8b)       :: totalflip,avclsize    

      real(r8b),external         :: rlfsr113
      external                      lfsrinit

      character datfile*14,avenfile*14,avmafile*14,avcofile*14
      datfile='heis000000.dat'
      avenfile='aven000000.dat'
      avmafile='avma000000.dat'
      avcofile='avco000000.dat'


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Start of main program
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      N_IMPSITE=nint(L*L*impconc)

      print *,'Program heiswolff5'
      print *,'--------------------------------------------------'
      print *,'L= ', L, ' Lt= ', Lt
      print *,'MC steps: ', NEQ, ' + ', NMESS

      write(datfile(5:7),'(I3.3)') L
      write(avenfile(5:7),'(I3.3)') L
      write(avmafile(5:7),'(I3.3)') L
      write(avcofile(5:7),'(I3.3)') L
      write(datfile(8:10),'(I3.3)') Lt
      write(avenfile(8:10),'(I3.3)') Lt
      write(avmafile(8:10),'(I3.3)') Lt
      write(avcofile(8:10),'(I3.3)') Lt

      open(2,file=datfile,status='replace')
      rewind(2)
      write(2,*) 'program heiswolff5'
      write(2,*) 'linear system sizes, L=', L, ' Lt= ', Lt
      write(2,*) 'equilibration steps     ',  NEQ      
      write(2,*) 'measurement steps   ', NMESS
      write(2,*) 'impurity conc.  ', real(impconc),',  number of imps. ',N_IMPSITE
      write(2,*) 'disorder configurations ', NCONF 
      write(2,*) 'CANON_DIS ',CANON_DIS
      write(2,*) 'LFSR-Init        ', IRINIT
      write(2,*) '-----------------'
      write(2,*) '  T     magnetization     energy   spec.heat  susceptibility   Binder  ncluster avclsize xit xis'
      close(2)


! Set up neigbor table

      do i1=0, Lt-1
      do i2=0, L-1
      do i3=0, L-1
          is = L*L*i1 + L*i2 + i3

          if (i1.eq.Lt-1) then
              m1(is)=is - L*L*(Lt-1)
          else
              m1(is)=is + L*L
          endif
          
          if (i1.eq.0) then
              m2(is)=is + L*L*(Lt-1)
          else
              m2(is)=is - L*L
          endif
          
          if (i2.eq.L-1) then
              m3(is)= is - L*(L-1)
          else
              m3(is)= is + L
          endif
          
          if (i2.eq.0) then
              m4(is)=is + L*(L-1)
          else
              m4(is)=is - L
          endif

          if (i3.eq.L-1) then
              m5(is)= is - (L-1)
          else
              m5(is)= is + 1
          endif

          if (i3.eq.0) then
              m6(is)= is + (L-1)
          else
              m6(is)= is - 1
          endif
      enddo
      enddo
      enddo


      confmag(:)   = 0.D0
      conf2mag(:)  = 0.D0
      confmag2(:)  = 0.D0
      confmag4(:)  = 0.D0
      conflogmag(:)= 0.D0
      confsusc(:)  = 0.D0
      confbin(:)    = 0.D0
      conf2susc(:) = 0.D0
      conf2bin(:)   = 0.D0
      confen(:)    = 0.D0
      confsph(:)   = 0.D0
      conf2en(:)   = 0.D0
      conf2sph(:)  = 0.D0
      confGt(:)    = 0.D0 
      confGs(:)    = 0.D0 
      confGtcon(:) = 0.D0 
      confGscon(:) = 0.D0 
      confxit(:)   = 0.D0
      confxis(:)   = 0.D0
      confxitcon(:)= 0.D0
      confxiscon(:)= 0.D0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Loop over disorder configurations   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      disorder: do iconf=1,NCONF 

      print *, 'disorder configuration ', iconf


! Initialize random number generator, IRINIT must be positive
      init=irinit+iconf-1
      call lfsrinit(init)

! Initialize impurities
      occu(:)=.true.
      iimp=0
      IF (CANON_DIS) THEN
      do while (iimp.lt.N_IMPSITE)
         ispace = int(rlfsr113()*L*L)
         if (occu(ispace)) then
            do itime=0,Lt-1
               occu(ispace+L*L*itime)=.false.
!               print *,ispace+L*L*itime 
            enddo
            iimp=iimp+1
         endif
      enddo
      ELSE
      do ispace=0,L*L-1
         if (rlfsr113()<IMPCONC) then
            do itime=0,Lt-1
               occu(ispace+L*L*itime)=.false.
            enddo
            iimp=iimp+1
         endif 
      enddo
      ENDIF 

      open(2,file=datfile,position='append') 
      write(2,*) 'disorder configuration ', iconf, ' imputities: ', iimp 
      close(2) 

! Initialize spins

      do is=0, L3-1
        if (occu(is)) then
          sx(is)= rlfsr113()-0.5D0
          sy(is)= rlfsr113()-0.5D0
          sz(is)= rlfsr113()-0.5D0
          slen=sqrt(sx(is)*sx(is) + sy(is)*sy(is) +sz(is)*sz(is))
          sx(is)=sx(is)/slen
          sy(is)=sy(is)/slen
          sz(is)=sz(is)/slen
        endif
      enddo

! Loop over temperatures
      ncluster=NCLUSTER0
      T=TMAX
      itemp=1 
      temperature: do while (T.gt.TMIN)
      beta=1./T
      print *,'T= ',real(T)

!     Equilibration, carry out NEQ full MC steps

      do isweep=1,NEQ/2
!         call metro_sweep
         call wolff_sweep(nsweepflip) 
      enddo
      totalflip=0     
      do isweep=1, NEQ/2
!         call metro_sweep
         call wolff_sweep(nsweepflip) 
         totalflip=totalflip+nsweepflip
      enddo
      avclsize=totalflip/ncluster/(NEQ-20)
      ncluster=int(L3/avclsize+0.5)+2
    
! Measuring loop, carry out NMESS full MC sweeps

      mag= 0.D0
      mag2=0.D0
      mag4=0.D0
      en=  0.D0
      magqt=0.D0
      magqs=0.D0
      Gt=  0.D0
      Gs=  0.D0
      
      totalflip=0
      do isweep=1, NMESS

!      call metro_sweep
      call wolff_sweep(nsweepflip)
      totalflip=totalflip+nsweepflip 

      call measurement     
      call corr_func

      enddo     ! of do isweep ...

      mag=mag/NMESS
      mag2=mag2/NMESS
      mag4=mag4/NMESS
      en=en/NMESS
      en2=en2/NMESS
      magqt=magqt/NMESS
      magqs=magqs/NMESS
      Gt=Gt/NMESS
      Gs=Gs/NMESS

      susc=(mag2 - mag**2)*L3*beta
      bin=1-mag4/(3*mag2**2)
      sph     =(en2-en**2)*L3*beta**2
      Gtcon=Gt-magqt**2
      Gscon=Gs-magqs**2

!      print *, mag2-mag**2, Gtcon

      xit= (mag2 - Gt)/ (Gt*qtime*qtime)
      xit=sqrt(abs(xit))
      xis= (mag2 - Gs)/ (Gs*qspace*qspace)
      xis=sqrt(abs(xis))
      xitcon= ((mag2 - mag**2) - Gtcon)/ (Gtcon*qtime*qtime)
      xitcon=sqrt(abs(xitcon))
      xiscon= ((mag2 - mag**2) - Gscon)/ (Gscon*qspace*qspace)
      xiscon=sqrt(abs(xiscon))

      avclsize=(totalflip/ncluster)/NMESS

      open(2,file=datfile,position='append') 
      write(2,'(1x,6(e12.6,1x),I6,10(1x,e12.6))') t,mag,en,sph,susc,bin,ncluster,avclsize, xit,xis
      close(2)

      confmag(itemp)=confmag(itemp)      +mag
      conf2mag(itemp)=conf2mag(itemp)    +mag**2
      confmag2(itemp)=confmag2(itemp)    +mag2
      confmag4(itemp)=confmag4(itemp)    +mag4
      conflogmag(itemp)=conflogmag(itemp)+log(mag)
      confsusc(itemp)=confsusc(itemp)    +susc
      confbin(itemp)=confbin(itemp)      +bin
      conf2susc(itemp)=conf2susc(itemp)  +susc**2
      conf2bin(itemp)=conf2bin(itemp)    +bin**2
      confen(itemp)=confen(itemp)        +en
      confsph(itemp)=confsph(itemp)      +sph
      conf2en(itemp)=conf2en(itemp)      +en**2
      conf2sph(itemp)=conf2sph(itemp)    +sph**2
      confGt(itemp)=confGt(itemp)        +Gt
      confGs(itemp)=confGs(itemp)        +Gs
      confGtcon(itemp)=confGtcon(itemp)  +Gtcon
      confGscon(itemp)=confGscon(itemp)  +Gscon
      confxit(itemp)=confxit(itemp)       +xit
      confxis(itemp)=confxis(itemp)       +xis
      confxitcon(itemp)=confxitcon(itemp) +xitcon
      confxiscon(itemp)=confxiscon(itemp) +xiscon
  

      ncluster=int(L3/avclsize+0.5)+2
      itemp=itemp+1
      T=T+DT
      enddo temperature

      open(7,file=avenfile,status='replace')
      rewind(7) 
      write(7,*) 'program heiswolff5'
      write(7,*) 'spatial + temporal system size', L, Lt
      write(7,*) 'equilibration steps     ',  NEQ
      write(7,*) 'measurement steps   ', NMESS
      write(7,*) 'impurity conc.  ', real(impconc),',  number of imps. ',N_IMPSITE
      write(7,*) 'disorder configurations ', NCONF 
      write(7,*) 'disorder configurations processed ',iconf 
      write(7,*) 'CANON_DIS ',CANON_DIS
      write(7,*) 'LFSR-Init        ', IRINIT
      write(7,*) '-----------------'
      write(7,*) '  T    [<energy>]  std.en    [<spec.heat>]    std.spec.heat'
      T=TMAX      
      itemp=1
      do while (T.gt.TMIN)
         write(7,'(1x,25(e12.6,1x))') t,confen(itemp)/iconf,sqrt(conf2en(itemp)/iconf-(confen(itemp)/iconf)**2),&
                                      confsph(itemp)/iconf, sqrt(conf2sph(itemp)/iconf-(confsph(itemp)/iconf)**2)
         T=T+DT
         itemp=itemp+1
      enddo
      close(7)
      open(7,file=avmafile,status='replace')
      rewind(7) 
      write(7,*) 'program heiswolff5'
      write(7,*) 'spatial + temporal system size', L, Lt
      write(7,*) 'equilibration steps     ',  NEQ
      write(7,*) 'measurement steps   ', NMESS
      write(7,*) 'impurity conc.  ', real(impconc),',  number of imps. ',N_IMPSITE
      write(7,*) 'disorder configurations ', NCONF 
      write(7,*) 'disorder configurations processed ',iconf 
      write(7,*) 'CANON_DIS ',CANON_DIS
      write(7,*) 'LFSR-Init        ', IRINIT
      write(7,*) '-----------------'
      write(7,*) ' T    [<mag>]  [<mag>^2]  std.mag  [<susc>] std.susc  [<Binder>]  std.bin [log <mag>]  Global.Binder'
      T=TMAX      
      itemp=1
      do while (T.gt.TMIN)
         write(7,'(1x,25(e12.6,1x))') t,confmag(itemp)/iconf,conf2mag(itemp)/iconf,&
              sqrt(conf2mag(itemp)/iconf-(confmag(itemp)/iconf)**2),confsusc(itemp)/iconf,&
              sqrt(conf2susc(itemp)/iconf-(confsusc(itemp)/iconf)**2), confbin(itemp)/iconf,&
              sqrt(conf2bin(itemp)/iconf-(confbin(itemp)/iconf)**2), conflogmag(itemp)/iconf,&
              1-confmag4(itemp)/iconf/(3*(confmag2(itemp)/iconf)**2)
         T=T+DT
         itemp=itemp+1
      enddo
      close(7)
      open(7,file=avcofile,status='replace')
      rewind(7) 
      write(7,*) 'program heiswolff5'
      write(7,*) 'spatial + temporal system size', L, Lt
      write(7,*) 'equilibration steps     ',  NEQ
      write(7,*) 'measurement steps   ', NMESS
      write(7,*) 'impurity conc.  ', real(impconc),',  number of imps. ',N_IMPSITE
      write(7,*) 'disorder configurations ', NCONF 
      write(7,*) 'disorder configurations processed ',iconf 
      write(7,*) 'CANON_DIS ',CANON_DIS
      write(7,*) 'LFSR-Init        ', IRINIT
      write(7,*) '-----------------'
      write(7,*) ' T  [xit]  [xis] [xit]/Lt  [xis]/L [xitcon]  [xiscon] [xitcon]/Lt  [xiscon]/L'
      write(7,*) ' glxit glxis  glxit/Lt glxis/L glxitcon glxiscon glxitcon/Lt glxiscon/L'
      T=TMAX      
      itemp=1
      do while (T.gt.TMIN)
         glxit= (confmag2(itemp) - confGt(itemp))/ (confGt(itemp)*qtime*qtime)
         glxit=sqrt(abs(glxit))
         glxis= (confmag2(itemp) - confGs(itemp))/ (confGs(itemp)*qspace*qspace)
         glxis=sqrt(abs(glxis))
         glxitcon= ((confmag2(itemp) - (confmag(itemp)**2)/iconf) - confGtcon(itemp))/ (confGtcon(itemp)*qtime*qtime)
         glxitcon=sqrt(abs(glxitcon))
         glxiscon= ((confmag2(itemp) - (confmag(itemp)**2)/iconf) - confGscon(itemp))/ (confGscon(itemp)*qspace*qspace)
         glxiscon=sqrt(abs(glxiscon))

         write(7,'(1x,25(e12.6,1x))') t,confxit(itemp)/iconf, confxis(itemp)/iconf,&
              confxit(itemp)/iconf/Lt, confxis(itemp)/iconf/L,&
              confxitcon(itemp)/iconf,confxiscon(itemp)/iconf,confxitcon(itemp)/iconf/Lt,confxiscon(itemp)/iconf/L,&
              glxit,glxis,glxit/Lt,glxis/L,glxitcon,glxiscon,glxitcon/Lt,glxiscon/L
         T=T+DT
         itemp=itemp+1
      enddo
      close(7)


      enddo disorder

      stop


      contains

      subroutine wolff_sweep(nflip)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!    Performs a Wolff sweep consisting of several single cluster flips
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      integer(i4b)    :: stack(0:L3-1)         ! stack for cluster construction
      integer(i4b)    :: sp                    ! stackpointer
      logical(ilog)   :: oldspin,newspin       ! temp. spin variables
      integer(i4b)    :: current               ! current site in cluster construction
      integer(i4b)    :: isize                 ! size of current cluster     
      integer(i4b)    :: nflip                 ! number of flipped spins
      integer(i4b)    :: icluster              ! cluster index
      real(r8b)       :: nx,ny,nz              ! reflection direction
      real(r8b)       :: scalar1,scalar2, snsn          ! scalar products in addition probability
      real(r8b)       :: padd,help
      
      nflip=0
      icluster=0
      
      do while(icluster.lt.NCLUSTER)
         is = int(L3*rlfsr113())             ! seed site for Wolff cluster
         if (occu(is)) then                  ! is site occupied?

         stack(0)=is
         sp=1

         nx= rlfsr113()-0.5D0
         ny= rlfsr113()-0.5D0
         nz= rlfsr113()-0.5D0
         slen= sqrt(nx*nx+ny*ny+nz*nz)
         nx= nx/slen
         ny= ny/slen
         nz= nz/slen
    
         icluster=icluster+1
         isize=1
         scalar2= (nx*sx(is)+ny*sy(is)+nz*sz(is))               ! scalar product for p_add

         sx(is)=sx(is)-2*nx*scalar2                               ! now flip seed spin
         sy(is)=sy(is)-2*ny*scalar2 
         sz(is)=sz(is)-2*nz*scalar2 

         do while(sp.gt.0)                   ! now build the cluster
           sp=sp-1
           current = stack(sp)              ! get site from stack

               scalar1=- (nx*sx(current)+ny*sy(current)+nz*sz(current))                ! scalar product for p_add

               if (occu(m6(current))) then
               scalar2=  (nx*sx(m6(current))+ny*sy(m6(current))+nz*sz(m6(current)) )   ! scalar product for p_add
               snsn=scalar1*scalar2
               if (snsn>0) then
                  padd=1.D0-exp(-2*beta*snsn)                                                 ! check whether parallel
                  help=rlfsr113()
                  if(help.lt.padd) then
                     sx(m6(current))=sx(m6(current))-2*nx*scalar2                               ! now flip current spin
                     sy(m6(current))=sy(m6(current))-2*ny*scalar2 
                     sz(m6(current))=sz(m6(current))-2*nz*scalar2 
                     stack(sp)=m6(current)
                     sp=sp+1
                     isize=isize+1
                  endif
               endif
               endif

               if (occu(m5(current))) then
               scalar2=  (nx*sx(m5(current))+ny*sy(m5(current))+nz*sz(m5(current)) )   ! scalar product for p_add
               snsn=scalar1*scalar2
               if (snsn>0) then
                  padd=1.D0-exp(-2*beta*snsn)                                                 ! check whether parallel
                  help=rlfsr113()
                  if(help.lt.padd) then
                     sx(m5(current))=sx(m5(current))-2*nx*scalar2                               ! now flip current spin
                     sy(m5(current))=sy(m5(current))-2*ny*scalar2 
                     sz(m5(current))=sz(m5(current))-2*nz*scalar2 
                     stack(sp)=m5(current)
                     sp=sp+1
                     isize=isize+1
                  endif
               endif
               endif

               if (occu(m4(current))) then
               scalar2=  (nx*sx(m4(current))+ny*sy(m4(current))+nz*sz(m4(current)) )   ! scalar product for p_add
               snsn=scalar1*scalar2
               if (snsn>0) then
                  padd=1.D0-exp(-2*beta*snsn)                                                 ! check whether parallel
                  help=rlfsr113()
                  if(help.lt.padd) then
                     sx(m4(current))=sx(m4(current))-2*nx*scalar2                               ! now flip current spin
                     sy(m4(current))=sy(m4(current))-2*ny*scalar2 
                     sz(m4(current))=sz(m4(current))-2*nz*scalar2 
                     stack(sp)=m4(current)
                     sp=sp+1
                     isize=isize+1
                  endif
               endif
               endif

               if (occu(m3(current))) then
               scalar2=  (nx*sx(m3(current))+ny*sy(m3(current))+nz*sz(m3(current)) )   ! scalar product for p_add
               snsn=scalar1*scalar2
               if (snsn>0) then
                  padd=1.D0-exp(-2*beta*snsn) 
                  help=rlfsr113()
                  if(help.lt.padd) then
                     sx(m3(current))=sx(m3(current))-2*nx*scalar2                               ! now flip current spin
                     sy(m3(current))=sy(m3(current))-2*ny*scalar2 
                     sz(m3(current))=sz(m3(current))-2*nz*scalar2 
                     stack(sp)=m3(current)
                     sp=sp+1
                     isize=isize+1
                  endif
               endif
               endif

                if (occu(m2(current))) then
               scalar2=  (nx*sx(m2(current))+ny*sy(m2(current))+nz*sz(m2(current)) )   ! scalar product for p_add
               snsn=scalar1*scalar2
               if (snsn>0) then
                  padd=1.D0-exp(-2*beta*snsn)                                                 ! check whether parallel
                  help=rlfsr113()
                  if(help.lt.padd) then
                     sx(m2(current))=sx(m2(current))-2*nx*scalar2                               ! now flip current spin
                     sy(m2(current))=sy(m2(current))-2*ny*scalar2 
                     sz(m2(current))=sz(m2(current))-2*nz*scalar2 
                     stack(sp)=m2(current)
                     sp=sp+1
                     isize=isize+1
                  endif
               endif
               endif

               if (occu(m1(current))) then
               scalar2=  (nx*sx(m1(current))+ny*sy(m1(current))+nz*sz(m1(current)) )   ! scalar product for p_add
               snsn=scalar1*scalar2
               if (snsn>0) then
                  padd=1.D0-exp(-2*beta*snsn)                                                 ! check whether parallel
                  help=rlfsr113()
                  if(help.lt.padd) then
                     sx(m1(current))=sx(m1(current))-2*nx*scalar2                               ! now flip current spin
                     sy(m1(current))=sy(m1(current))-2*ny*scalar2 
                     sz(m1(current))=sz(m1(current))-2*nz*scalar2 
                     stack(sp)=m1(current)
                     sp=sp+1
                     isize=isize+1
                  endif
               endif
               endif

            enddo       ! of cluster building and flipping

         nflip   = nflip   +isize

         endif          ! of if(occu(is))
      enddo          ! of do icluster

      end subroutine wolff_sweep


     subroutine metro_sweep
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    carries out one Metropolis sweep
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      do is=0, L3-1
      if (occu(is)) then
          nsx=0.D0
          nsy=0.D0
          nsz=0.D0
          if (occu(m1(is))) then 
             nsx=nsx+sx(m1(is))
             nsy=nsy+sy(m1(is))
             nsz=nsz+sz(m1(is))
          endif
          if (occu(m2(is))) then 
             nsx=nsx+sx(m2(is))
             nsy=nsy+sy(m2(is))
             nsz=nsz+sz(m2(is))
          endif
          if (occu(m3(is))) then 
             nsx=nsx+sx(m3(is))
             nsy=nsy+sy(m3(is))
             nsz=nsz+sz(m3(is))
          endif
          if (occu(m4(is))) then 
             nsx=nsx+sx(m4(is))
             nsy=nsy+sy(m4(is))
             nsz=nsz+sz(m4(is))
          endif
          if (occu(m5(is))) then 
             nsx=nsx+sx(m5(is))
             nsy=nsy+sy(m5(is))
             nsz=nsz+sz(m5(is))
          endif
          if (occu(m6(is))) then 
             nsx=nsx+sx(m6(is))
             nsy=nsy+sy(m6(is))
             nsz=nsz+sz(m6(is))
          endif

          sxnew= rlfsr113()-0.5D0
          synew= rlfsr113()-0.5D0
          sznew= rlfsr113()-0.5D0
          slen=sqrt(sxnew*sxnew + synew*synew +sznew*sznew)
          sxnew=sxnew/slen
          synew=synew/slen
          sznew=sznew/slen

          dE= nsx*(sx(is)-sxnew)+ nsy*(sy(is)-synew)+nsz*(sz(is)-sznew)
          if (dE<0.or.(exp(-dE*beta)>rlfsr113())) then
             sx(is)=sxnew
             sy(is)=synew
             sz(is)=sznew
          endif
      endif     ! of if(occu(is))
      enddo     ! of do is ...


      end subroutine metro_sweep


      subroutine measurement
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! measures energy and magnetization
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      sweepen=0
      mx=0
      my=0
      mz=0
      do is=0,L3-1
         if (occu(is)) then
            mx=mx+sx(is)
            my=my+sy(is)
            mz=mz+sz(is)
            nsx=0
            nsy=0
            nsz=0
            if (occu(m1(is))) then
               nsx=nsx+sx(m1(is))
               nsy=nsy+sy(m1(is))
               nsz=nsz+sz(m1(is))
            endif
            if (occu(m3(is))) then
               nsx=nsx+sx(m3(is))
               nsy=nsy+sy(m3(is))
               nsz=nsz+sz(m3(is))
            endif
            if (occu(m5(is))) then
               nsx=nsx+sx(m5(is))
               nsy=nsy+sy(m5(is))
               nsz=nsz+sz(m5(is))
            endif
            sweepen=sweepen-sx(is)*nsx
            sweepen=sweepen-sy(is)*nsy
            sweepen=sweepen-sz(is)*nsz
         endif
      enddo
      en=en+sweepen/L3
      en2=en2+(sweepen/L3)**2
      sweepmag=   sqrt(mx*mx+my*my+mz*mz)/L3
      mag= mag +  sweepmag 
      mag2=mag2+  sweepmag**2 
      mag4=mag4+  sweepmag**4

      end subroutine measurement

      subroutine corr_func
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! measures correlation function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      real(r8b)   :: Remqtx,Remqty,Remqtz,Immqtx,Immqty,Immqtz
      real(r8b)   :: Remq2x,Remq2y,Remq2z,Immq2x,Immq2y,Immq2z
      real(r8b)   :: Remq3x,Remq3y,Remq3z,Immq3x,Immq3y,Immq3z
 
      Remqtx=0
      Remqty=0
      Remqtz=0
      Immqtx=0
      Immqty=0
      Immqtz=0
     
      Remq2x=0
      Remq2y=0
      Remq2z=0
      Immq2x=0
      Immq2y=0
      Immq2z=0

      Remq3x=0
      Remq3y=0
      Remq3z=0
      Immq3x=0
      Immq3y=0
      Immq3z=0

      do i1=0, Lt-1
      do i2=0, L-1
      do i3=0, L-1
         is = L*L*i1 + L*i2 + i3
         if (occu(is)) then
            Remqtx=Remqtx+sx(is)*cos(qtime*i1)
            Remqty=Remqty+sy(is)*cos(qtime*i1)
            Remqtz=Remqtz+sz(is)*cos(qtime*i1)
            Immqtx=Immqtx+sx(is)*sin(qtime*i1)
            Immqty=Immqty+sy(is)*sin(qtime*i1)
            Immqtz=Immqtz+sz(is)*sin(qtime*i1)

            Remq2x=Remq2x+sx(is)*cos(qspace*i2)
            Remq2y=Remq2y+sy(is)*cos(qspace*i2)
            Remq2z=Remq2z+sz(is)*cos(qspace*i2)
            Immq2x=Immq2x+sx(is)*sin(qspace*i2)
            Immq2y=Immq2y+sy(is)*sin(qspace*i2)
            Immq2z=Immq2z+sz(is)*sin(qspace*i2)

            Remq3x=Remq3x+sx(is)*cos(qspace*i3)
            Remq3y=Remq3y+sy(is)*cos(qspace*i3)
            Remq3z=Remq3z+sz(is)*cos(qspace*i3)
            Immq3x=Immq3x+sx(is)*sin(qspace*i3)
            Immq3y=Immq3y+sy(is)*sin(qspace*i3)
            Immq3z=Immq3z+sz(is)*sin(qspace*i3)
         endif
      enddo
      enddo
      enddo
      sweepmagqt=Remqtx**2+Remqty**2+Remqtz**2 + Immqtx**2+Immqty**2+Immqtz**2
      sweepmagqt=sqrt(sweepmagqt)/L3
      magqt=magqt+sweepmagqt
      Gt=Gt+sweepmagqt**2

      sweepmagq2=(Remq2x**2+Remq2y**2+Remq2z**2 + Immq2x**2+Immq2y**2+Immq2z**2)
      sweepmagq3=(Remq3x**2+Remq3y**2+Remq3z**2 + Immq3x**2+Immq3y**2+Immq3z**2)
      sweepmagq2=sqrt(sweepmagq2)/L3
      sweepmagq3=sqrt(sweepmagq3)/L3
      magqs=magqs+0.5*(sweepmagq2+sweepmagq3)
      Gs=Gs+0.5*(sweepmagq2**2+sweepmagq3**2)

      end subroutine corr_func
      end


!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
! Random number generator RLFSR113  - real version
!
! Following a suggestion of Pierre L'Ecuyer 1997
! "Tables of maximally equidistributed combined LFSR generators"
! see http://www.iro.umontreal.ca/~lecuyer/papers.html
!
! A call to rlfsr113() gives one random real in the open
! interval (0,1).
!
! Before using rlfsr113 call lfsrinit(seed) to initialize
! the generator by random integers produced by Park/Millers
! minimal standard LCG.
! Seed should be any positive integer.
! 
! FORTRAN version by Thomas Vojta, vojta@physik.tu-chemnitz.de
! 
! History:
!  04 Feb 1998    v0.9    first FORTRAN implementation
!  05 Feb 1998    v0.91   corrected multiplicator am to 1/(2^31)
!  15 Apr 1999    v0.92   added real*8 rlfs113 in lfsrinit
!  10 Feb 2000    v0.93   changed lfsrinit to allow for CONSTANT arguments
! 
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      FUNCTION rlfsr113 ()
      implicit none
      integer,parameter         :: i4b= kind(2147483647)
      integer,parameter         :: r8b= kind(1.D200)

      real(r8b),parameter       :: am=4.656612873077d-10

      real(r8b)    rlfsr113             
      integer(i4b) b,z1,z2,z3,z4
      common /lfsrcom/z1,z2,z3,z4

      b  = ishft(ieor(ishft(z1,6),z1),-13)
      z1 = ieor(ishft(iand(z1,-2),18),b)

      b  = ishft(ieor(ishft(z2,2),z2),-27)
      z2 = ieor(ishft(iand(z2,-8),2),b)

      b  = ishft(ieor(ishft(z3,13),z3),-21)
      z3 = ieor(ishft(iand(z3,-16),7),b)

      b  = ishft(ieor(ishft(z4,3),z4),-12)
      z4 = ieor(ishft(iand(z4,-128),13),b)

      rlfsr113=ishft( ieor(ieor(ieor(z1,z2),z3),z4) , -1)*am
!      print *, 'RLFSR113 ', rlfsr113

      return
      end


      SUBROUTINE lfsrinit(iinit)
      implicit none
      integer,parameter         :: r8b= kind(1.D200)
      integer,parameter         :: i4b= kind(2147483647)

      integer(i4b) idum,ia,im,iq,ir,iinit
      integer(i4b) k,z1,z2,z3,z4,c1,c2,c3,c4
      real(r8b)    rlfsr113,rdum
      parameter (ia=16807,im=2147483647,iq=127773,ir=2836)
      common /lfsrcom/z1,z2,z3,z4

      data c1 /B'11111111111111111111111111111110'/
      data c2 /B'11111111111111111111111111111000'/
      data c3 /B'11111111111111111111111111110000'/
      data c4 /B'11111111111111111111111110000000'/
      if ((c1.ne.-2).or.(c2.ne.-8).or.(c3.ne.-16).or.(c4.ne.-128)) then
         print *,'Nonstandard integer representation. Stoped.'
         stop
      endif

      idum=iinit
      if (idum.le.0) idum=1
      k=(idum)/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum = idum + IM
      if (idum.lt.2) then
         z1=idum+2 
      else 
         z1=idum
      endif
      k=(idum)/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum = idum + IM
      if (idum.lt.8) then 
         z2=idum+8 
      else 
         z2=idum
      endif
      k=(idum)/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum = idum + IM
      if (idum.lt.16) then
         z3=idum+16 
      else 
         z3=idum
      endif
      k=(idum)/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum = idum + IM
      if (idum.lt.128) then
         z4=idum+128 
      else 
         z4=idum
      endif

      rdum=rlfsr113()
      
      return
      end