! FORTRAN 90 VERSION  
!       (Compile  with  
!             f90 -free anneal.f        DECs
!             xlf90 anneal.f            IBM
!             f90 +source=free anneal.f HP  (f90 anneal.f90 also works)
!             f90  -freeform anneal.f   SGI 
!             f90 -free anneal.f        SUN  )
PROGRAM anneal       
  IMPLICIT NONE
  INTEGER :: nsteps=1000                    !number of steps in a block
  INTEGER, PARAMETER :: ndim=1
  INTEGER, PARAMETER :: ENOUT=11, XOUT=12
  REAL*8, DIMENSION(ndim)  :: x, xsum, x2sum, xmin
  REAL*8  :: fx
  REAL*8  :: ensum, en2sum, enav, en2av, cv !variables to compute specific heat
  REAL*8  :: temperature
  REAL*8  :: tfactor                        ! Temperature Reduction Factor
  INTEGER :: naccepted
  REAL*8  :: aratio, anorm, emin, sigmax
  INTEGER :: istep
  REAL*8, EXTERNAL ::  func
  INTEGER :: rand_seed(18), seed_size !!Variables for setting the seed

  !  Simulated annealing program for minimization of a fct.
  !  using straightforward Metropolis algorithm.
  !
  ! Provide: 
  !   use GNUPLOT graphical plot to elaborate results.
  !   details of Annealing Schedule and iterations to success.
  !   use "specific heat" to determine if annealing schedule should be slowed.


  !Open files for energy and postion record.
  OPEN(UNIT=ENOUT,FILE="en.out",STATUS="UNKNOWN")
  OPEN(UNIT=XOUT,FILE="x.out",STATUS="UNKNOWN")

  !Set starting point for searching for minimum.
  x = 1.0; fx=func(x); emin=fx
  WRITE(6,'("#Starting from x=",5f10.5)') x
  WRITE(6,'("#with function value f=",f10.5)') fx


  !Define your annealing schedule in this loop.
! temperature = 10d0
  WRITE(6,*) '  input initial (high) temperature (e.g., 10)'
  READ(5,*) temperature
  if(temperature.lt. 0.001) then
      temperature=0.001
     WRITE(6,*) '  temperature reset =',temperature
  endif
! tfactor = 0.9d0                              !Factor for annealing schedule.
  WRITE(6,*) '  input annealing factor (e.g., 0.9)'
  READ(5,*) tfactor
  WRITE(6,*) '  number of steps per block (equilibriation, e.g., 1000)'
  READ(5,*) nsteps
  !Set the random seed.
  CALL RANDOM_SEED(SIZE=seed_size)
  CALL RANDOM_SEED(GET=rand_seed(1:seed_size))
  WRITE(*,'("Enter ",i2," integer seed for random &
           &number generator, current seed is",18i14)')  &
           seed_size,rand_seed(1:seed_size)
  READ(*,*) rand_seed(1:seed_size)
  CALL RANDOM_SEED(PUT=rand_seed(1:seed_size))

  !Before starting loop, thermalize for block
  DO istep=1,nsteps
    CALL step(x(1:1),fx,temperature,func,naccepted)
  END DO
  WRITE(*,'("#Themalized for ",i5," steps.")') nsteps


  WRITE(ENOUT,'(a8,2x,a11,3(2x,a8),a11)') &
       "#temp ","<E>  ", "C_v  ", "sigma_x  ", "ac.ratio","E_min"
  
  WRITE(XOUT,'("#temp, then ndim cols for <x>, and ndim cols for x_min")')
  
  anneal_cycle: DO

    !First a block metropolis steps
    naccepted = 0 
    !We keep track of <E> and <E^2> to compute C_v
    ensum=0d0; en2sum=0d0; xsum=0d0; x2sum=0d0
    DO istep=1,nsteps
      CALL step(x(1:ndim),fx,temperature,func,naccepted)
      ensum = ensum+fx
      en2sum = en2sum+fx**2
      xsum = xsum + x
      x2sum = x2sum + x**2
      IF (fx<emin) THEN
        emin = fx
        xmin = x
      END IF
    END DO

    anorm=1d0/FLOAT(nsteps)
    enav = ensum*anorm
    cv = (en2sum*anorm-enav**2)/temperature**2
    sigmax = SQRT(SUM( x2sum*anorm - (xsum*anorm)**2 ))
    aratio=naccepted*anorm

    WRITE(ENOUT,'(e8.2,2x,e11.5,3(2x,e8.2),2x,e11.5)') &
       temperature, enav, cv, sigmax, aratio, emin 

    WRITE(XOUT,'(11(e12.6,2x))') &
       temperature, x(1:ndim), xmin(1:ndim)

    !Anneal by decreasing the temperature
    temperature = temperature * tfactor

    IF (temperature < 10d-5) EXIT anneal_cycle
  END DO anneal_cycle

  CLOSE(XOUT)
  CLOSE(ENOUT)
END PROGRAM anneal

SUBROUTINE step(x,fx,temperature,funct,naccepted)
  !Perform a Metropolis step.
  IMPLICIT NONE
  !Argument variables
  REAL*8, INTENT(INOUT), DIMENSION(1) :: x
  REAL*8, INTENT(INOUT) :: fx
  REAL*8, INTENT(IN) :: temperature
  REAL*8, EXTERNAL :: funct
  INTEGER, INTENT(INOUT) :: naccepted
  !Local variables
  REAL*8 :: x_new(SIZE(x))
  REAL*8 :: fx_new
  REAL*8, EXTERNAL :: rand
  REAL*8, PARAMETER :: scale=0.5      ! SCALE should be chosen for specific FUNC

  x_new  = x+ scale*SQRT(temperature)*(rand() - 0.5) ! stochastic moves in range
  fx_new = funct(x_new)                               ! new object function value
  IF ( EXP(-(fx_new-fx)/temperature) > rand() ) THEN
    naccepted = naccepted+1 
    fx = fx_new               ! if successful, save
    x(:) = x_new(:)
  END IF 
END SUBROUTINE step


! Random Number Generator  
REAL*8 FUNCTION rand() 
  REAL*8 a
  CALL RANDOM_NUMBER(a)
  rand=a
END FUNCTION rand


! Funtion to Minimize
REAL*8 FUNCTION func(x)
  IMPLICIT NONE
  REAL*8, INTENT(IN), DIMENSION(1) :: x
  func = (  (x(1) + 0.2 ) * x(1) + cos( 14.5 * x(1) - 0.3 ) ) 
!       *(  (y + 0.5 ) * y + cos( 14.5 * y - 0.4 ) )
END FUNCTION
