/* (Quantenmechanische) MC-Simulation eines 'Quantendots'
   Numerische Methoden der Vielteilchenphysik WS01/02
   Roland Winkler und Guenter Zwicknagel */

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

#define PI 3.14159265359

int np, anzmon;
double omega2, coul, temp;

double u_eff (double **partx, double **party);
double upot (double *partx, double *party);

void force (double *x, double *y, double **f);
double dforce ();

double ham_op (double **partx, double **party, double **f);

float ran1 (int *idum);
float gasdev (int *idum);

void mswap (double ***x, double ***y);
double **dmatrix(int nrl, int nrh, int ncl, int nch);

/*  #define ZUFALL (ran1 (&idum) - 0.5) */
#define ZUFALL gasdev (&idum)

/*****************/
/* Hauptprogramm */
/*****************/

int main (int argc, char *argv[])
{
  extern int np, anzmon;
  extern double omega2, coul, temp;
  int anzges, idum, i, j, k, schritt, akzep;
  double zeta, rskala, upot_alt, upot_neu, du, ham, ham_sum;
  double **partx, **party, **xn, **yn, **frc;
  int startsamp = 20000, anzmess = 1000;
  FILE *daten_file;

  if (argc == 7) {
    sscanf (argv [1], "%lf", &omega2);   /* omega2            */
    sscanf (argv [2], "%lf", &coul);     /* coul              */
    sscanf (argv [3], "%lf", &temp);     /* Temperatur        */
    sscanf (argv [4], "%i",  &np);       /* Teilchenzahl      */
    sscanf (argv [5], "%i",  &anzmon);   /* Zahl der Monomere */
    sscanf (argv [6], "%i",  &anzges);   /* Zahl der Schritte */
  }
  else {
    printf ("%s: Wrong number of arguments\n", argv [0]);
    exit (1);
  }

  /* Felder alloziieren */
  partx = dmatrix (0, anzmon, 1, np);
  party = dmatrix (0, anzmon, 1, np);
  xn    = dmatrix (0, anzmon, 1, np);
  yn    = dmatrix (0, anzmon, 1, np);

  frc   = dmatrix (1, np, 1, 2); /* Hilsfeld fuer Erwartungswert <H> */

  /* Die charakteristische Laengenskala */
  rskala = 4.0 / sqrt (omega2) + 5.0 * pow (coul/omega2, 1.0/3.0);

  /* Random-Generator initialisieren */
  idum = - fabs (rskala*np*temp + 123.) - 0.5;
  printf ("Initialisierung: %f\n", ran1 (&idum));

  /* Teilchenorte initialisieren */
  zeta = 0.3 * rskala;
  for (j = 1; j <= np; j++) {
    partx [0][j] = zeta * ZUFALL;
    party [0][j] = zeta * ZUFALL;
  }
  /* Idee: Die einzelnen Monomere liegen im Phasenraum eng zusammen */
  zeta = 0.01 * zeta;
  for (i = 1; i <= anzmon; i++) {
    for (j = 1; j <= np; j++) {
      partx [i][j] = partx [0][j] + zeta * ZUFALL;
      party [i][j] = party [0][j] + zeta * ZUFALL;
    }
  }

  /* Berechnung des Potentials fuer die Anfangskonfiguration */
  upot_alt = u_eff (partx, party);
  ham = ham_op (partx, party, frc);
  printf ("Startpotential   %e\n", ham);

  /* Ausgabe initialisieren */
  daten_file = fopen ("datqdclass.dat", "w");
  fprintf (daten_file, "# omega2:            %f\n", omega2);
  fprintf (daten_file, "# coul:              %f\n", coul);
  fprintf (daten_file, "# temp:              %f\n", temp);
  fprintf (daten_file, "# Teilchenzahl:      %i\n", np);
  fprintf (daten_file, "# Zahl der Monomere: %i\n", anzmon);
  fprintf (daten_file, "# Zahl der Schritte: %i\n", anzges);

  akzep = 0;
  zeta  = 0.02 * rskala;
  ham_sum = 0.0;

  k = 0;

  /* Beginn des Schleifenprogramms */
  for (schritt = 1; schritt <= anzges; schritt++) {

    /* Erzeugen einer neuen Konfiguration               */
    /* Idee: Variiere abwechselnd immer nur ein Monomer */
    for (i = 0; i <= anzmon; i++) {
      for (j = 1; j <= np; j++) {
        xn [i][j] = partx [i][j];
        yn [i][j] = party [i][j];
      }
    }
    k = (k+1) % (anzmon+1);
    for (j = 1; j <= np; j++) {
      xn [k][j] = partx [k][j] + zeta * ZUFALL;
      yn [k][j] = party [k][j] + zeta * ZUFALL;
    }

    /* Berechnung der Aenderung der pot. Energie */
    upot_neu = u_eff (xn, yn);
    du = upot_neu - upot_alt;

    /* Befrage das Orakel */
    if (exp (-du / temp) > ran1 (&idum)) {
      akzep ++;
      upot_alt = upot_neu;
      mswap (&partx, &xn);
      mswap (&party, &yn);
      /* Nur wenn sich die Konfiguration geaendert hat, wird H ausgewertet */
      ham = ham_op (partx, party, frc);
    }

    /* Observable bestimmen */
    if (schritt > startsamp) {
      ham_sum += ham;
    }

    if ((schritt % anzmess) == 0) {
      /* Datenausgabe */
      if (schritt > startsamp) {
        fprintf (daten_file, "%i %13.6e  %13.6e  %13.6f\n",
                 schritt,
                 ham_sum / (np * (schritt - startsamp)),
                 (double) akzep / anzmess, zeta);
      }

      /* Anpassung des zeta-Parameters */
      if ((double) akzep / anzmess > 0.5)
        zeta *= 1.05;
      else
        zeta *= 0.95;

      akzep = 0;
    }
  }

  printf ("Letztes Potential   %e\n", ham);

  return 0;
}

/************************/
/* Effektives Potential */
/************************/

double u_eff (double **partx, double **party)
{
  extern int np, anzmon;
  extern double temp;
  double u, upoly, rx, ry, efak, pfak;
  int i, j;

  /* Einmal haben wir immer das echte Potential */
  pfak = 1.0 / (anzmon+1);
  u    = pfak * upot (partx [0], party[0]);

  if (anzmon > 0) {
    efak = 0.5 * (anzmon+1) * temp * temp;
    for (i = 1; i <= anzmon; i++) {
      /* Harmonische Polymer-Wechselwirkung */
      upoly = 0.0;
      for (j = 1; j <= np; j++) {
        rx = partx [i][j] - partx [i-1][j];
        ry = party [i][j] - party [i-1][j];
        upoly += rx * rx + ry * ry;
      }
      u += efak * upoly + pfak * upot (partx [i], party[i]);
    }

    /* Ring-Polymer */
    upoly = 0.0;
    for (j = 1; j <= np; j++) {
      rx = partx [0][j] - partx [anzmon][j];
      ry = party [0][j] - party [anzmon][j];
      upoly += rx * rx + ry * ry;
    }
    u += efak * upoly;
  }

  return u;
}

/************************/
/* Echte Wechselwirkung */
/************************/

double upot (double *partx, double *party)
{
  extern int np;
  extern double omega2, coul;
  double u = 0.;
  int i, j;

  /* harmonischer Anteil */
  for (i = 1; i <= np; i++)
    u += partx [i] * partx [i] + party [i] * party [i];
  u *= 0.5 * omega2;

  /* Coulomb-Anteil */
  if (coul != 0.0) {
    for (i = 1; i <= np-1; i++) {
      for (j = i+1; j <= np; j++) {
        double rx = partx [i] - partx [j];
        double ry = party [i] - party [j];
        u += coul / sqrt (rx*rx + ry*ry);
      }
    }
  }
  return u;
}

/*********************/
/* Hamilton-Operator */
/*********************/

double ham_op (double **partx, double **party, double **frc)
{
  extern int np, anzmon;
  extern double temp;
  double fak = (anzmon+1) * temp;
  double ham = 0.0;
  int i;

  if (anzmon > 0) {
    /* Kraefte aus Polymer-Wechselwirkung */
    for (i = 1; i <= np; i++) {
      frc [i][1] = fak * (partx [0][i] - partx [1][i]);
      frc [i][2] = fak * (party [0][i] - party [1][i]);
    }

    /* Kraefte der echten Wechselwirkung draufaddieren */
    force (partx [0], party [0], frc);

    /* Quadrieren und aufsummieren */
    for (i = 1; i <= np; i++)
      ham += frc [i][1] * frc [i][1] + frc [i][2] * frc [i][2];
  }

  ham = pow ((anzmon+1) * temp / (4 * PI * PI), (double) np * (anzmon+1))
    * (- 0.5 * (ham + np * ( - 2.0 * fak - dforce () / fak))
       + upot (partx [0], party [0]));

  return ham;
}

/***********/
/* Kraefte */
/***********/

void force (double *x, double *y, double **frc)
{
  extern int np, anzmon;
  extern double omega2, coul;
  int i, j;
  double fx, fy, rx, ry, sq, den;
  double ofak = omega2 * (anzmon+1) * temp;
  double cfak = coul / ((anzmon+1) * temp);

  /* Kraefte beim harmonischen Oszillator */
  for (i = 1; i <= np; i++) {
    frc [i][1] += ofak * x[i];
    frc [i][2] += ofak * y[i];
  }

  /* Kraefte aus Coulomb-Wechselwirkung */
  if (coul != 0.0) {
    for (i = 1; i <= np; i++) {
      for (j = 1; j < i; j++) {
        rx  = x [i] - x [j];
        ry  = y [i] - y [j];
        sq  = sqrt (rx * rx + ry * ry);
        den = cfak / (sq * sq * sq);
        fx  = rx * den;
        fy  = ry * den;
        /* Actio = Reactio */
        frc [i][1] -= fx;
        frc [i][2] -= fy;
        frc [j][1] += fx;
        frc [j][2] += fy;
      }
    }
  }
  return;
}

/***********************************/
/* Zweite Ableitung des Potentials */
/***********************************/

double dforce ()
{
  extern double omega2;
  return 2.0 * omega2;
}

/*****************/
/* Swap matrices */
/*****************/

void mswap (double ***x, double ***y)
{
  double **temp;
  temp = *x;
  *x   = *y;
  *y   = temp;
}

/***********************************************************************/
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
/***********************************************************************/

#define NR_END 1

double **dmatrix(int nrl, int nrh, int ncl, int nch)
{
  int i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  double **m;

  /* allocate pointers to rows */
  m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
  m += NR_END;
  m -= nrl;

  /* allocate rows and set pointers to them */
  m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
  m[nrl] += NR_END;
  m[nrl] -= ncl;

  for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

  /* return pointer to array of pointers to rows */
  return m;
}

#undef NR_END

/**************************************************/
/* Zufallszahlengenerator aus "Numercial Recipes" */
/**************************************************/

#define M1 259200
#define IA1 7141
#define IC1 54773
#define RM1 (1.0/M1)
#define M2 134456
#define IA2 8121
#define IC2 28411
#define RM2 (1.0/M2)
#define M3 243000
#define IA3 4561
#define IC3 51349

float ran1(int *idum)
{
  static long ix1,ix2,ix3;
  static float r[98];
  float temp;
  static int iff=0;
  int j;

  if (*idum < 0 || iff == 0) {
    iff=1;
    ix1=(IC1-(*idum)) % M1;
    ix1=(IA1*ix1+IC1) % M1;
    ix2=ix1 % M2;
    ix1=(IA1*ix1+IC1) % M1;
    ix3=ix1 % M3;
    for (j=1;j<=97;j++) {
      ix1=(IA1*ix1+IC1) % M1;
      ix2=(IA2*ix2+IC2) % M2;
      r[j]=(ix1+ix2*RM2)*RM1;
    }
    *idum=1;
  }
  ix1=(IA1*ix1+IC1) % M1;
  ix2=(IA2*ix2+IC2) % M2;
  ix3=(IA3*ix3+IC3) % M3;
  j=1 + ((97*ix3)/M3);
  temp=r[j];
  r[j]=(ix1+ix2*RM2)*RM1;
  return temp;
}

#undef M1
#undef IA1
#undef IC1
#undef RM1
#undef M2
#undef IA2
#undef IC2
#undef RM2
#undef M3
#undef IA3
#undef IC3

/*********************************/
/* Gauss-verteilte Zufallszahlen */
/*********************************/

float gasdev (int *idum)
{
  static int iset=0;
  static float gset;
  float fac,r,v1,v2;

  if  (iset == 0) {
    do {
      v1=2.0*ran1(idum)-1.0;
      v2=2.0*ran1(idum)-1.0;
      r=v1*v1+v2*v2;
    } while (r >= 1.0);
    fac=sqrt(-2.0*log(r)/r);
    gset=v1*fac;
    iset=1;
    return v2*fac;
  } else {
    iset=0;
    return gset;
  }
}