/***********************************************************/
/* PROGRAM:    qmc.c                                       */
/*                                                         */
/* AUTHOR      Adriaan Schakel                             */
/* DATE:       January, 2007                               */
/* SYNOPSIS:   Feynman Path Integral Monte Carlo           */
/* REFERENCE:  P. K. MacKeown, Evaluation of Feynman       */
/*             path integrals by Monte Carlo methods,      */
/*             Am. J. Phys. 53, 880 (1985)                 */
/* COMMENTS:   Adapted from R.H. Landau and M.J. Paez,     */
/*             Computational Physics (Wiley-Interscience,  */
/*             New York, N.Y., 1997)                       */
/***********************************************************/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

      
#define MCSMAX 1000000		/* maximum number of MC sweeps        */
#define NEQ    200000		/* number of equilibration sweeps     */
#define SEED   68111
#define DELTA  1.5              /* maximum trial change               */
#define DTAU   1.               /* delta tau                          */  
#define N      400              /* number of "atoms" = Trotter number */
#define M      400              /* number of bins                     */
#define XMAX   4.               /* maximum value of x                 */  

/* values of DELTA and DTAU should be chosen so that about 50% of
   updates are accepted */

int
main()
{
   double en, x[N], xo, xmax;
   int i, j, k, n, r, psi2[M], mcs;
   
   double energy (double xo, double xk, double xkp1, double xkm1);
                /* gives change in energy for update */
   
   FILE *output;                           /* save data in qmc.dat */
   output = fopen("qmc_08.dat","w");
   
   srand48(time(0));
 /*   srand48(SEED); */			   /* seed number generator */

   xmax = 0.;
   mcs = 0;
   for (j=0; j<N; j++) 
     x[j]=0.;     /* initial path */
   

   for (j=0; j<M; j++) 
     psi2[j]=0;    /* and psi2 */

   r = 0;
   for (i=0; i<MCSMAX; i++) {
     k = drand48()*N;             /* choose an atom at random */

     xo = x[k];
     x[k] += (2.*drand48()-1.)*DELTA;   /* trial displacement */


     /* energy change */

     en = energy(xo,x[k],x[(k+1)%N],x[(k-1+N)%N]);

     /* % is to ensure periodic boundary conditions  */
     
     if ((en>0) && (exp(-DTAU*en) < drand48())) {
       x[k] = xo;	   /* reject */
       r++;               /* number of rejections */
     }
 
     
     if (i > NEQ) {  /* start measuring after equilibration sweeps */ 
       mcs++;
       n = M*x[k]/(2*XMAX) + M/2;
       psi2[n]++;
     }
   }

   for (n=0; n<M; n++)  {
     fprintf(output, "%f\t%f\n", (double) 2.*n*XMAX/M - XMAX, 
	     (double) M/2/XMAX*psi2[n]/mcs);
     /*      fprintf(output, "%d\t%f\n", n, x[n]);  */
   }
   printf("data stored in qmc_08.dat\n");
   printf("rejection rate: %f\n",(float) r/MCSMAX);
   fclose(output);

   return (0);
}


/* function returns change in energy for update */

double energy (double xo, double xk, double xkp1, double xkm1)
{
  double delta;
  
  delta= ((xkp1-xk)*(xkp1-xk) + (xkm1-xk)*(xkm1-xk) -
          (xkp1-xo)*(xkp1-xo) - (xkm1-xo)*(xkm1-xo))/(2.*DTAU*DTAU) 
          + (xk*xk - xo*xo)/2.;   /* any other potential can be taken */

   return (delta);
}