c ***************************************************************************
c Copyright(C) M.Kotelyanskii, November, 1997
c This code is just an example, it can be freely distributed and modified
c as soon as this notice remains in the file.
c It has been tested on the IBM Risc6000 under AIX. Other systems may require
c changes, particularly the functions calls for the random numbers generation.
c Even though it works properly to the best of the author's knowledge,
c he does not assume any responsibility for any damage caused by it
c ***************************************************************************
c ***************************************************************************
c *****  This is example of the Monte Carlo method for a 
c ***** non--self-avoiding polymer chain on the lattice                              
c ***** To make things interesting we look at the polar chain in the electric
c ***** field along the positive X direction. 
c ***** Each bond has a dipole moment one.
c ***************************************************************************

       implicit real*8(a-h,o-z)

c **** Specify parameters:
c **** Electric field strength
       parameter(Field = 8.0)
c **** Chain length
       parameter(Length=100)
c **** Number of trials
       parameter(Ncycle=10000, Nprint=Ncycle/10)
c **** Temperature
       parameter(T=0.1, beta=1./T)
c **** Random number generator's seed
       parameter(iseed=1234567)

c **** set counters
c **** Here we will calculate and average the end-to-end vector 
       rxav  = 0.0d0
       rx2av = 0.0d0
       ryav  = 0.0d0
       ry2av  = 0.0d0
       ree2av  = 0.0d0

c *** and energy

      Eav=0
      E2av=0

c *** Precompute the weights for going along x,-x,y,-y
      wpx = exp( beta*Field)
      wnx = exp(-beta*Field)
      wpy = 1.0 
      wny = 1.0 

      sw=(wpx+wnx+wpy+wny)

      wpx=wpx/sw               
      wpy=wpy/sw               
      wnx=wnx/sw               
      wny=wny/sw               

c *** Place them at the interval from 0 to 1 
c *** for decision making

      wnx=wnx+wpx
      wpy=wnx+wpy
      wny=wpy+wny

c *** set seed for the random numbers generator (notice, that the name srand
c *** may be different on different architectures )
       call srand(iseed)

c *** Start MC cycle (assuming rand() returns random real*1 from 0 to 1)

       do icycle = 1, Ncycle

c *** We always start from the origin 
         x=0.0
         y=0.0
         Energy = 0.0
c *** Loop through chain segments building the chain
        do isegment=1,Length
c ***   get the next segment's position, and energy, associated with it
c ***   Next segment can go back, forward, left, or right.
c **** decide which way to build the next bond 
          r=rand()
          if(r.lt.wpx) then
c **** we go to the +x
             x=x+1.0              
             Energy = Energy - E
          else if(r.lt.wnx) then
c **** we go to the -x
             x=x-1.0              
             Energy = Energy + E
             else if(r.lt.wpy) then
c **** we go to +y, no energy change
                y=y+1.0              
                else if(r.lt.wny) then
c **** we go to -y, no energy change
                y=y-1.0              
           endif

c **** end building chain loop
         enddo

c **** increment averages
c ****  average the end-to-end vector 
       rxav    = rxav   + x 
       rx2av   = rx2av  + x*x 
       ryav    = ryav   + y 
       ry2av   = ry2av  + y*y 
       ree2av  = ree2av + x*x+y*y 

c *** and energy

      Eav=Eav+Energy
      E2av=E2av+Energy*Energy

c **** print out the running averages
           if(mod(icycle,Nprint).eq.1) then
             rxrunav = rxav / icycle
             rx2runav = rx2av / icycle
             rx2runav = sqrt((rx2runav-rxrunav*rxrunav)/icycle)
             ryrunav = ryav / icycle
             ry2runav = ry2av / icycle
             ry2runav = sqrt((ry2runav-ryrunav*ryrunav)/icycle)
             ree2runav = ree2av / icycle
             Erunav=Eav/icycle
             E2runav=E2runav/icycle
             E2runav = sqrt((E2runav-Erunav*Erunav)/icycle)
             write(6,2) icycle,Erunav,E2runav,rxrunav,rx2runav,ryrunav,
     >         ry2runav,ree2runav 
2      format(1x,i8, "cycles Average Energy = ",e12.5," Std.Dev ",e12.5,
     >        1x, "Rx = ",e12.5, " Std.Dev ", e12.5,
     >        1x, "Ry = ",e12.5, " Std.Dev ", e12.5,
     >        1x, "Ree2 = ",e12.5)
           endif

c *** End the cycles loop
       enddo

c *** That's all folks

        end

