/* cc -O2 -o alg42 alg42.c -lm

  Simulace Lennardovy-Jonesovy tekutiny v periodickych okrajovych
  podminkach s vypoctem energie metodou nejblizsiho obrazu (viz
  poznamka na str. 57 skript).
  Program podle opravene verze algoritmu 4.2, str. 47 skript.

 */

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

#define N 200 /* pocet atomu */
double seed=1; /* nasada pro generator nahodnych cisel */
const double d=0.4; /* polovina delky zkusebniho posunuti */
const double T=1.2; /* teplota (v redukovanych jednotkach) */
const double rho=0.8; /* hustota (v redukovanych jednotkach) */
const double k=1; /* Boltzmannova konstanta v redukovanych jednotkach */

typedef struct { double x,y,z; } vector;

double rnd(void)
/* nahodne cislo v intervalu [0,1) */
{
  seed=fmod(seed*16807e0,2147483647e0);
  return seed/2147483647e0;
}

double L; /* delka hrany simulacni bunky */
double Lh; /* L/2 */

double u(vector r1,vector r2) 
/* Lennarduv-Jonesuv potencial v redukovanych jednotkach */
{
  vector dr;
  double rr;

  /* viz algoritmus 5.2: r1 i r2 musi byt v zakladni bunce */
  dr.x=Lh-fabs(Lh-fabs(r1.x-r2.x));
  dr.y=Lh-fabs(Lh-fabs(r1.y-r2.y));
  dr.z=Lh-fabs(Lh-fabs(r1.z-r2.z));

  rr=dr.x*dr.x+dr.y*dr.y+dr.z*dr.z; /* kvadrat vzdalenosti nejblizsich sousedu */
  rr*=rr*rr;

  return (4/rr-4)/rr;
}

int main(int narg,char **arg)
{
  int i,j;
  vector r[N]; /* konfigurace */
  double Epair[N][N]; /* tabulka parovych energii */
  double Ezkus[N]; /* tabulka parovych energii zkusebni konfigurace */
  double Utot,Uzkus,U;

  L=pow(N/rho,1./3); /* hrana simulacni bunky */
  Lh=L/2;

  /* pocatecni konfigurace je nahodna */
  for (i=0; i<N; i++) {
    r[i].x=rnd()*L;
    r[i].y=rnd()*L;
    r[i].z=rnd()*L;
  }

  /* naplneni tabulky parovych energii */
  for (i=0; i<N-1; i++)
    for (j=i+1; j<N; j++) /* nyni (i,j) prochazi vsechny pary castic, i<j */
      Epair[i][j]=Epair[j][i]=u(r[i],r[j]);
  
  for (;;) { /* forever... */

    int pocet_prijatych=0;

    Utot=0;

    for (i=0; i<N; i++) {
      vector rzkus;

      /* zkusebni posunuti atomu i */
      rzkus.x=r[i].x+d*(rnd()-0.5);
      rzkus.y=r[i].y+d*(rnd()-0.5);
      rzkus.z=r[i].z+d*(rnd()-0.5);

      /* a normalizace do zakladni bunky */
      if (rzkus.x<0) rzkus.x+=L; else if (rzkus.x>L) rzkus.x-=L;
      if (rzkus.y<0) rzkus.y+=L; else if (rzkus.y>L) rzkus.y-=L;
      if (rzkus.z<0) rzkus.z+=L; else if (rzkus.z>L) rzkus.z-=L;
      
      /* vypocet stare (U) a nove (Uzkus) energie */
      Uzkus=U=0;
      /* soucet pres vsechny atomy ruzne od i */
      for (j=0; j<N; j++) if (i!=j) {
	Uzkus+=Ezkus[j]=u(r[j],rzkus); /* nove parove energie schovame */
	U+=Epair[i][j];
      }

      /* Metropolisuv test */
      if (rnd()<exp(-(Uzkus-U)/(k*T))) { /* konfigurace prijata */
	pocet_prijatych++;
	r[i]=rzkus;
	Utot+=Uzkus;
	for (j=0; j<N; j++) /* aktualizace tabulky parovych energii */
	  Epair[i][j]=Epair[j][i]=Ezkus[j];
      }

      else /* konfigurace odmitnuta */
	Utot+=U;

    }
    printf("chi=%4.2f U/atom=%5.3f\n",
	   (double)pocet_prijatych/N, Utot/(2*N));
  }
  return 0;
}