      program metrop_2N_1D

c     Metropolis calculation
c     for 2N system in 1D
c     with initial state equal to 1 
c     for all spatial lattice sites

      implicit integer(i-n)
      implicit double precision(a-h,o-y)
      implicit complex*16(z)

cccccc

      parameter(myseed = 510271)
      parameter(numtrials = 100000000)
      parameter(nwarmup = 100000)
      parameter(nprintevery = 100000)
      parameter(Lt = 50)
      parameter(L = 6)
      parameter(ainv = 100.D0)
      parameter(atinv = 100.D0)
      parameter(alphat = ainv/atinv)
      parameter(c = -0.2D0)
      parameter(am = 938.92D0)
      dimension v(0:L-1,0:Lt)
      dimension dualv(0:L-1,0:Lt)
      dimension s(0:L-1,0:Lt-1)
      dimension temp(0:L-1)
      dimension overlap(0:Lt)

cccccc

      call sgrnd(myseed)

      accept = 0.D0
      ratio_bin = 0.D0

c     since we have the same initial vector for both up
c     and down spins and the auxiliary-field transfer matrix 
c     is independent of spin, we can use the same single-particle
c     states for up and down spins

      do nx = 0,L-1
         v(nx,0) = 1.D0/dsqrt(1.D0*L)
         dualv(nx,Lt) = 1.D0/dsqrt(1.D0*L)
      enddo     

      bose = 0.D0
      do nt = 0,Lt-1
         do nx = 0,L-1
            s(nx,nt) = grnd()-0.5D0
            bose = bose + 0.5D0*s(nx,nt)*s(nx,nt)
         enddo
      enddo

      call getv(v,s,0,Lt,c,am,alphat,L,Lt)
      amp = 0.D0
      do nx = 0,L-1
         amp = amp + dualv(nx,Lt)*v(nx,Lt)
      enddo
     
      do ntrial = 1,numtrials

         call getdualv(dualv,s,0,Lt,c,am,alphat,L,Lt)
         
         do nt = 0,Lt-1
            do nx = 0,L-1
               s_old = s(nx,nt)
               s_new = s(nx,nt) + grnd() - 0.5D0            
               bosediff = 0.5D0*(s_new*s_new-s_old*s_old)
               ampnew = amp      
     $              + dualv(nx,nt+1)*v(nx,nt)
     $              *dsqrt(-c*alphat)*(s_new-s_old)
               if (grnd() .lt. 
     $              (ampnew/amp)**2.D0*dexp(-bosediff)) then
                  accept = accept + 1.D0               
                  amp = ampnew
                  s(nx,nt) = s_new
                  bose = bose + bosediff
               endif
            enddo
            call getv(v,s,nt,nt+1,c,am,alphat,L,Lt)
         enddo

         mtrial = ntrial - nwarmup
         
         if (mtrial .gt. 0) then
            amp1 = 0.D0
            do nx = 0,L-1
               amp1 = amp1 + dualv(nx,Lt)*v(nx,Lt-1)
            enddo            
            ratio_bin = ratio_bin + amp1/amp
            energy = -dlog(mtrial/ratio_bin)/alphat         
            if (mod(mtrial,nprintevery) .eq. 0) then
               write(*,*)
               write(*,*)'mtrial',mtrial,'energy (MeV)',energy*ainv
               write(*,*)'accept',accept/(ntrial*L*Lt)
            endif            
         endif

      enddo         
         
      stop
      end


      subroutine getv(v,s,ntm,ntp,c,am,alphat,L,Lt)

      implicit integer(i-n)
      implicit double precision(a-h,o-y)
      implicit complex*16(z)
      dimension v(0:L-1,0:Lt)      
      dimension s(0:L-1,0:Lt-1)      

      do nt = ntm+1,ntp
         do nx = 0,L-1
            v(nx,nt) = 
     $           v(nx,nt-1)*(1.D0 - 1.D0/am*alphat
     $           + dsqrt(-c*alphat)*s(nx,nt-1))
     $           + v(mod(nx-1+L,L),nt-1)/(2.D0*am)*alphat
     $           + v(mod(nx+1,L),nt-1)/(2.D0*am)*alphat
         enddo
      enddo

      return
      end


      subroutine getdualv(dualv,s,ntm,ntp,c,am,alphat,L,Lt)

      implicit integer(i-n)
      implicit double precision(a-h,o-y)
      implicit complex*16(z)
      dimension dualv(0:L-1,0:Lt)      
      dimension s(0:L-1,0:Lt-1)      

      do nt = ntp,ntm+1,-1
         do nx = 0,L-1
            dualv(nx,nt-1) = 
     $           dualv(nx,nt)*(1.D0 - 1.D0/am*alphat
     $           + dsqrt(-c*alphat)*s(nx,nt-1))
     $           + dualv(mod(nx-1+L,L),nt)/(2.D0*am)*alphat
     $           + dualv(mod(nx+1,L),nt)/(2.D0*am)*alphat
         enddo
      enddo

      return
      end

! ************************************************************************

! A C-program for MT19937: Real number version
!   genrand() generates one pseudorandom real number (double)
! which is uniformly distributed on [0,1]-interval, for each
! call. sgenrand(seed) set initial values to the working area
! of 624 words. Before genrand(), sgenrand(seed) must be
! called once. (seed is any 32-bit integer except for 0).
! Integer generator is obtained by modifying two lines.
!   Coded by Takuji Nishimura, considering the suggestions by
! Topher Cooper and Marc Rieffel in July-Aug. 1997.
!
! This library is free software; you can redistribute it and/or
! modify it under the terms of the GNU Library General Public
! License as published by the Free Software Foundation; either
! version 2 of the License, or (at your option) any later
! version.
! This library is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
! See the GNU Library General Public License for more details.
! You should have received a copy of the GNU Library General
! Public License along with this library; if not, write to the
! Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
! 02111-1307  USA
!
! Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
! When you use this, send an email to: matumoto@math.keio.ac.jp
! with an appropriate reference to your work.
!
! ************************************************************************
! Fortran translation by Hiroshi Takano.  Jan. 13, 1999.
!
!   genrand()      -> double precision function grnd()
!   sgenrand(seed) -> subroutine sgrnd(seed)
!                     integer seed
!
! This program uses the following non-standard intrinsics.
!   ishft(i,n): If n>0, shifts bits in i by n positions to left.
!               If n<0, shifts bits in i by n positions to right.
!   iand (i,j): Performs logical AND on corresponding bits of i and j.
!   ior  (i,j): Performs inclusive OR on corresponding bits of i and j.
!   ieor (i,j): Performs exclusive OR on corresponding bits of i and j.
!
! ************************************************************************

      subroutine sgrnd(seed)
      implicit integer(a-z)

! Period parameters
      parameter(N     =  624)
      dimension mt(0:N-1)
! the array for the state vector
      common /block/mti,mt
      save   /block/

!      setting initial seeds to mt[N] using
!      the generator Line 25 of Table 1 in
!      [KNUTH 1981, The Art of Computer Programming
!         Vol. 2 (2nd Ed.), pp102]

      mt(0)= iand(seed,-1)
      do 1000 mti=1,N-1
        mt(mti) = iand(69069 * mt(mti-1),-1)
 1000 continue

      return
      end

! ************************************************************************

      double precision function grnd()
      implicit integer(a-z)

! Period parameters
      parameter(N     =  624)
      parameter(N1    =  N+1)
      parameter(M     =  397)
      parameter(MATA  = -1727483681)
! constant vector a
      parameter(UMASK = -2147483648)
! most significant w-r bits
      parameter(LMASK =  2147483647)
! least significant r bits
! Tempering parameters
      parameter(TMASKB= -1658038656)
      parameter(TMASKC= -272236544)

      dimension mt(0:N-1)
! the array for the state vector
      common /block/mti,mt
      save   /block/
      data   mti/N1/
! mti==N+1 means mt[N] is not initialized

      dimension mag01(0:1)
      data mag01/0, MATA/
      save mag01

! mag01(x) = x * MATA for x=0,1

      TSHFTU(y)=ishft(y,-11)
      TSHFTS(y)=ishft(y,7)
      TSHFTT(y)=ishft(y,15)
      TSHFTL(y)=ishft(y,-18)

      if(mti.ge.N) then
! generate N words at one time
        if(mti.eq.N+1) then
! if sgrnd() has not been called,
          call sgrnd(4357)
! a default initial seed is used
        endif

        do 1000 kk=0,N-M-1
            y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
            mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1)))
 1000   continue
        do 1100 kk=N-M,N-2
            y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
            mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1)))
 1100   continue
        y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK))
        mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1)))
        mti = 0
      endif

      y=mt(mti)
      mti=mti+1
      y=ieor(y,TSHFTU(y))
      y=ieor(y,iand(TSHFTS(y),TMASKB))
      y=ieor(y,iand(TSHFTT(y),TMASKC))
      y=ieor(y,TSHFTL(y))

      if(y.lt.0) then
        grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0)
      else
        grnd=dble(y)/(2.0d0**32-1.0d0)
      endif

      return
      end

! ************************************************************************

