(0) Place WP1.tar.gz into a directory and

 $ tar -zxvf WP1.tar.gz
  
 Command will generate directory WP1 with subdirectories

  /source
  /runnable

(1) Compiling the 1D examples

 $ cd WP1
 $ cd source
 $ make clean
 $ make
 $ cd ../

(2) Running the 1D examples

 $ cd runnable
   (you may modify the "input" file to perform desired simulations)
 $ ./mcmc.x input > output

   Each run generates the following files:

       esurf: energy surface E(L) we would like to explore.
              structure of file: [ L  E ]

       adist: the analytical distribution(s), Pa(k) ~ exp(-E(L)/T[k]) at temperature
              T(k), T(k) is the temperature of the kth replica and k = 1,...,K, K being
              the total number of replicas.
              structure of file: [ L Pa(1) ... Pa(K) ]

       ndist: the corresponding numerical distributions, Pn(k) obtained by running the
              simulation where the input parameters are provided by students
              structure of the file: [ L Pn(1) ... Pa(K) ]

       energy: the energy of replica(s), E(k) as a function of the MC iteration steps
              structure of the file: [ i E(1) ... E(K) ]

       convd: distance metric, d(k) = |Pa(k) - Pn(k)|, measuring how much each replica's
              numerical distribution deviates from the analytical distribution as a function
              of the MC iteration steps
              structure of the file: [ i d(1) ... d(K) ]

(3) Description of input parameters

  --------------------------------------------------------------------------------------
  1: ~sim_gen_def[
  2:    \method{PT}            PT: Parallel Tempering Monte Carlo, EEMC: Equi-Energy Monte Carlo
  3:    \energy_surface{DWP}   DWP: Double Well Potential, RFP: Fourier Series Potential
  4:    \replica_number{1}     N >=1 N: Natural number ! N = 1 MC: Standard Monte Carlo
  5:    \total_step_mc{100000} N >=1
  6:    \local_step_md{10}     N >=1
  7:    \time_step_md{0.02}    R > 0   R: Real number
  8:    \statistics_freq{1000} N >=1
  9:    \prob_eemc_jump{0.15}  R in (0,1)
  10:   \temperature{1.0}      R > 0 Usually in (0.4,1.5)
  11:   \energy_gap{1.1}       R in (1,2) Usually in (1.1,1.2)
  12:   \random_seed{-2039480938l}
  13: ]
  ---------------------------------------------------------------------------------------

 - The three algorithms you could try: a) Monte Carlo, \method{PT} and \replica_number{1}
                                       b) Parallel Tempering Monte Carlo, \method{PT} and \replica_number{>1}
                                       c) Equi energy Monte Carlo, \method{EEMC} (*)
                                         (*) we recommend that you mainly focus on comparing a)
                                             and b) and perhaps run a few simulations with c)

 - Two energy surfaces provided:       a) Double Well Potential \energy_surface{DWP}
                                       b) Random Fourier Series Potential \energy_surface{RFP}

 - Monte Carlo move size / acceptance rate
   The size of an MC steps is defined by (\local_step_md{})*(\time_step_md{}) as a gradient
   based few step length dynamical path to generate the next proposed Monte Carlo move.
   Before each MC move, the velocities are resampled according to the Maxwell-Boltzmann distribution.
   Increasing the MC step size, leads to larger moves but lower acceptance rate.

 - (Local) MC move size / inter replica exchange
   The local MC step vs. Multi-Canonical MC "jump": At each MC iteration, a probabilistic
   approach is used to either propagate a particular replica independently at its given
   temperature or attempt a replica exchange with an adjacent higher temperature replica.
   The probability to attempt a replica exchange is controlled by (\prob_eemc_jump{})
   where eemc stands for extended-ensemble monte carlo. Thus, the probability of independent
   replica propagation is (1.0 - \prob_eemc_jump{}).

 - Separation of temperatures / inter replica jumps
   The tempearature separation can be increased by altering the \energy_gap{}, usually in the
   range 1.1-1.2. Historically, first we generate initial energy levels using
   Ei(k) = (\energy_gap{})^(k), k = 1,...,K+1 next we use a linear transformation to map
   these energy levels to an final energy interval with fixed E(1) and E(K+1) so that the
   temperature levels are defined as T(k) = T0*(E(k+1)-E(k))/(E(2)-E(1)), k = 1,...,K, T0
   being target temperature (\temperature{}) of the first replica.

(4) Suggested exercises

   Given the target distribution of replica k, (*) Pa(k) ~ exp(-E(L)/T[k]), the rate/success
   of barrier crossing mainly depends on the ratio (E(L)/T[k]), therefore we only let you change
   T[k] (by changing \temperature) in order to increase the difficulty level of overcoming
   barriers. Once you set up a particular system, and some difficulty levels (e.g. by using a
   range of temperatures from 1.3 to 0.4 from easier to more difficult) you could explore the
   parameter and algorithmic space to gain experience how to run a more efficient simulation
   that results in a good match between the analytical and numerical distribution within a
   smaller number of steps. For example, you could monitor the convergence rate as a function
   of the number MC steps and make notes how it depends on the input parameters. Some suggested
   questions:

   a) What is the most efficient acceptance rate/associated step size? (use \prob_eemc_jump{0})
   b) What is the most efficient temperature separation/inter replica jump rate?
   c) How does the increase in the replica number help overcoming the energy barriers?
   d) Observe, visualize and interpret the quantitative / qualitative measures monitoring
      the success of exploring the model space.
   e) Make comparisons between the different energy surfaces studies and the algorithms
      used to approximate the known limiting (analytical distributions)
   d) You could also make distribution of different energy levels and plot them. Is there
      an overlap between adjacent energy distributions? How does it depend on the
      temperature separation.

                                               (*) Pa(k) ~ exp{-E(L)/(kb T[k])}, we use kb = 1.
                                                   kb : Boltzmann constant



    Hope you will learn and benefit from this practical exercise and we look forward to
    your questions, notes and interacting with you during the practical.