/****************************************************************/
/*** Simulates system of Lennard-Jones particles              ***/
/*** A.K. Hartmann May 2004                                   ***/
/*** Rechnergestütztes Wissenschaftliches Rechnen             ***/
/*** University of Goettingen, Germany 2004                   ***/
/****************************************************************/

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

/* stores data of one atom: */
typedef struct
{
  double        m;       /* mass of atom */
  double    sigma;       /* 'size' of atom for LJ */
  double  epsilon;       /* LJ energy parameter */
  double       *x;       /* position of atom */
  double       *v;       /* velocity of atom */
  double       *f;       /* force on atom */
} glas_atom_t;

/* stores all global system data: */
typedef struct
{
  int               dim;       /* dimension of system */
  double             *l;       /* sizes of system */
  double            *lh;       /* half sizes of system */
  int                 N;       /* total number of particles */
  double            tau;       /* integration step size */
  double           tau2;       /* integration step size^2 */
  double              T;       /* temperature */
  double         virial;       /* Virial = sum_(i>j) r_ij*f_ij/3 */ 
} glas_system_t;


/******************** glas_setup() **********************/
/** Creates and initialises atom data                  **/
/** PARAMETERS: (*)= return-paramter                   **/
/**     system: global system parameters               **/
/** RETURNS:                                           **/
/**     array of atom data                             **/ 
/********************************************************/
glas_atom_t *glas_setup(glas_system_t *system)
{
  glas_atom_t *atom;          /* here all atom data is stored */
  int t, t2, d;                              /* loop counters */

  atom = (glas_atom_t *) malloc(system->N*sizeof(glas_atom_t));
  for(t=0; t<system->N; t++)              /* initialize atoms */
  {
    atom[t].m = 1;
    atom[t].sigma = 1;
    atom[t].epsilon = 1;
    atom[t].x = (double *) malloc(system->dim*sizeof(double));
    atom[t].v = (double *) malloc(system->dim*sizeof(double));
    atom[t].f = (double *) malloc(system->dim*sizeof(double));
    for(d=0; d<system->dim; d++)         /* put atom randomly */
    {
      atom[t].x[d] = drand48()*system->l[d];
      atom[t].v[d] = 0;
    }
  }
  return(atom);
}


/***************** glas_setup_regular() *****************/
/** Puts atoms in regular way (almost cubic lattice)   **/
/** into system                                        **/
/** PARAMETERS: (*)= return-paramter                   **/
/**     system: global system parameters               **/
/**       atom: atom data                              **/ 
/**          n: how many atoms in each direction       **/
/** RETURNS:                                           **/
/**      nothing                                       **/
/********************************************************/
double glas_put_regular(glas_system_t *system, glas_atom_t *atom, int *n)
{
  int t, d;                                  /* loop counters */
  double dist,                  /* distance between two atoms */ 
    volume;
  int n0;                /* number of atoms in each direction */
  int x,y,z;                                /* coordinates */
  int dim;
  int N;                           /* number of atoms (local) */

  dim = system->dim;
  if(dim != 3)
  {
    fprintf(stderr, "glas_put_regular works only for 3d!\n");
    exit(1);
  }
  N = 1;
  for(d=0; d<dim; d++)
  {
    N*= n[d];
  }
  if(N != system->N)
  {
    fprintf(stderr, "Number of total atoms =%d != #atoms passed =%d\n",
	    system->N, N);
    exit(1);
  }
  t=0;
  for(z=0; z<n[2]; z++)
    for(y=0; y<n[1]; y++)
      for(x=0; x<n[0]; x++)
      {
	atom[t].x[0] = (x+0.5*(y%2)-0.5*(z%2))*system->l[0]/n[0];
	atom[t].x[1] = (y+0.33333*(z%2))*system->l[1]/n[1];
	atom[t].x[2] = z*system->l[2]/n[2];
	t++;
      }
}

/***************** glas_energy_forces() *****************/
/** Calculates potential energy and forces of system   **/
/** PARAMETERS: (*)= return-paramter                   **/
/**     system: global system parameters               **/
/**       atom: atom data                              **/ 
/** RETURNS:                                           **/
/**      energy                                        **/
/********************************************************/
double glas_energy_forces(glas_system_t *system, glas_atom_t *atom)
{
  double energy = 0.0;
  double epsilon4, sigma2;             /* interaction parameters */
  int t, t2, d;                                  /* loop counter */
  double *r;              /* difference vector between two atoms */
  double r2, rm6;                     /* distance^2, distance^-6 */
  double force;
  
  r = (double *) malloc(system->dim*sizeof(double));
  for(t=0; t<system->N; t++)             /* loop over all  atoms */
    for(d=0; d<system->dim; d++)
      atom[t].f[d] = 0;                      /* initialise force */
  system->virial = 0;

  for(t=0; t<system->N; t++)     /* loop over all pairs of atoms */
    for(t2=t+1; t2<system->N; t2++)
    {
      r2 = 0.0;                            /* calculate distance */
      d = 0;
      while(d<system->dim)
      {
	r[d] = atom[t].x[d] - atom[t2].x[d];
	if(r[d] < -system->lh[d])    /* minimum image convention */
	  r[d] += floor(-r[d]/system->l[d]+0.5)*system->l[d];
	if(r[d] > system->lh[d])
	  r[d] -= floor(r[d]/system->l[d]+0.5)*system->l[d];	  
	r2 += r[d]*r[d];
	d++;
      }
      sigma2 = 0.5*(atom[t].sigma+atom[t2].sigma);  
      sigma2 = sigma2*sigma2;
      epsilon4 = 4.0*sqrt(atom[t].epsilon * atom[t2].epsilon);
      r2 /= sigma2;
      rm6 = 1.0/(r2*r2*r2);
      energy += epsilon4*(rm6*(rm6-1.0));     /* calulate energy */
      for(d=0; d<system->dim; d++)
      {
	force = 6*epsilon4*(rm6*(2*rm6-1.0))*r[d]/(r2*sigma2);
	atom[t].f[d] += force;                     /* and forces */
	atom[t2].f[d] -= force;
	system->virial += r[d]*force;                  /* virial */
      }      
    }

  system->virial /= 3;                                 /* virial */
  free(r);
  return(energy);
}

/***************** glas_kinetic_energy() ****************/
/** Calculates kinetic energy of system                **/
/** PARAMETERS: (*)= return-paramter                   **/
/**     system: global system parameters               **/
/**       atom: atom data                              **/ 
/** RETURNS:                                           **/
/**      energy                                        **/
/********************************************************/
double glas_kinetic_energy(glas_system_t *system, glas_atom_t *atom)
{
  double energy = 0.0;
  int t, d;                                      /* loop counter */

  for(t=0; t<system->N; t++)              /* loop over all atoms */
    for(d=0; d<system->dim; d++)
      energy += 0.5*atom[t].m*atom[t].v[d]*atom[t].v[d];

  return(energy);
}

/***************** glas_rescale_kinetic() ***************/
/** Rescales velocity to yield temperature T           **/
/** PARAMETERS: (*)= return-paramter                   **/
/**     system: global system parameters               **/
/**       atom: atom data                              **/ 
/** RETURNS:                                           **/
/**      nothing                                       **/
/********************************************************/
void glas_rescale_kinetic(glas_system_t *system, glas_atom_t *atom)
{
  double energy = 0.0;
  double lambda;
  int t, d;                                      /* loop counter */

  energy = glas_kinetic_energy(system, atom)/system->N;
  lambda = sqrt(3.0*system->T/(2*energy));
  for(t=0; t<system->N; t++)              /* loop over all atoms */
    for(d=0; d<system->dim; d++)
      atom[t].v[d] *= lambda;                /* rescale velocity */
}


/***************** glas_velocity verlet() ***************/
/** Calculate new positions, velocities and forces     **/
/** using the velocity verlet algorithm                **/
/** Condition: current forces already calculated       **/ 
/** PARAMETERS: (*)= return-paramter                   **/
/**     system: global system parameters               **/
/**       atom: atom data                              **/ 
/** RETURNS:                                           **/
/**      potential energy                              **/
/********************************************************/
double glas_velocity_verlet(glas_system_t *system, glas_atom_t *atom)
{
  int t, d;                                      /* loop counter */
  double energy;                             /* potential energy */
  for(t=0; t<system->N; t++)              /* loop over all atoms */
    for(d=0; d<system->dim; d++)
    {
      atom[t].x[d] += atom[t].v[d]*system->tau  /* new positions */
	+0.5*atom[t].f[d]*system->tau2/atom[t].m;
      atom[t].v[d] += 0.5*atom[t].f[d]*system->tau/atom[t].m;
    }
  energy =  glas_energy_forces(system, atom);      /* new forces */
  for(t=0; t<system->N; t++)              /* loop over all atoms */
    for(d=0; d<system->dim; d++)
    {                                   /* finish new velocities */
      atom[t].v[d] += 0.5*atom[t].f[d]*system->tau/atom[t].m;
    }
  return(energy);

}


/***************** glas_plot_cfg() **********************/
/** Prints configuration 'atom'  and system boundaries **/
/** to povray file with name cfg.pov.                  **/
/** Display with x-povray +I cfg.pov                   **/
/** PARAMETERS: (*)= return-paramter                   **/
/**     system: global system parameters               **/
/**       atom: atom data                              **/ 
/** RETURNS:                                           **/
/**      nothing                                       **/
/********************************************************/
void glas_plot_cfg(glas_system_t *system, glas_atom_t *atom)
{
  char filename[1000];
  FILE *povfile;
  int t, d;                                   /* loop counters */
  double *r;                                       /* position */

  r = (double *) malloc(system->dim*sizeof(double));
  sprintf(filename, "cfg.pov");
  povfile = fopen(filename, "w");                 /* open file */
  fprintf(povfile, "#include \"colors.inc\"\n"   /* header etc */
	  "#include \"shapes.inc\"\n\n");
  fprintf(povfile, "background { color Yellow }\n\n");
  fprintf(povfile, "camera {\n  location <%f, %f, %f>\n", 
	  0.5*system->l[0], -1.5*system->l[1], 0.5*system->l[2]);
  fprintf(povfile, "  sky <0,0,1>\n");
  fprintf(povfile, "  look_at  <%f, %f, %f>\n}\n\n", 
	  0.5*system->l[0], 0.5*system->l[1], 0.5*system->l[2]);
  fprintf(povfile, 
	  "  light_source { <%f, %f, %f> color White}\n\n",
	  0.5*system->l[0], -0.5*system->l[1], 1.5*system->l[2]);
  fprintf(povfile, 
	  "  light_source { <%f, %f, %f> color White}\n\n",
	  -0.5*system->l[0], -0.5*system->l[1], 1.5*system->l[2]);

  fprintf(povfile,                   /* print system boundaries */
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  0.0, 0.0, 0.0, 
	  system->l[0], 0.0, 0.0);
  fprintf(povfile,       
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  system->l[0], 0.0, 0.0, 
	  system->l[0], system->l[1], 0.0);
  fprintf(povfile,       
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  system->l[0], system->l[1], 0.0, 
	  0.0, system->l[1], 0.0);
  fprintf(povfile,       
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  0.0, system->l[1], 0.0, 
	  0.0, 0.0, 0.0);
  fprintf(povfile,                  
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  0.0, 0.0, system->l[2], 
	  system->l[0], 0.0, system->l[2]);
  fprintf(povfile,       
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  system->l[0], 0.0, system->l[2], 
	  system->l[0], system->l[1], system->l[2]);
  fprintf(povfile,       
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  system->l[0], system->l[1], system->l[2], 
	  0.0, system->l[1], system->l[2]);
  fprintf(povfile,       
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  0.0, system->l[1], system->l[2], 
	  0.0, 0.0, system->l[2]);
  fprintf(povfile,       
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  0.0, 0.0, 0.0, 
	  0.0, 0.0, system->l[2]);
  fprintf(povfile,       
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  system->l[0], 0.0, 0.0, 
	  system->l[0], 0.0, system->l[2]);
  fprintf(povfile,       
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  system->l[0], system->l[1], 0.0, 
	  system->l[0], system->l[1], system->l[2]);
  fprintf(povfile,       
	  "cylinder{ <%f, %f, %f>, <%f, %f, %f>, 0.1\n"
	  "   pigment { Red } }\n",
	  0.0, system->l[1], 0.0, 
	  0.0, system->l[1], system->l[2]);


  for(t=0; t<system->N; t++)                     /* print atoms */
  {
    for(d=0; d<system->dim; d++)
    {
      r[d] = atom[t].x[d];
      if(r[d] < 0)                   /* fold positions into box */
	  r[d] += floor(-r[d]/system->l[d]+1)*system->l[d];
	if(r[d] > system->l[d])
	  r[d] -= floor(r[d]/system->l[d])*system->l[d];	  
    }
    if(system->dim == 3)
      fprintf(povfile, 
	      "sphere { <%f,%f,%f>, %f\n    pigment { Blue }}\n",
	      r[0], r[1], r[2], 0.4*atom[t].sigma);
    else if(system->dim == 2)
      fprintf(povfile, 
	      "sphere { <%f,%f,%f>, %f\n    pigment { Blue }}\n",
	      r[0], r[1], 0, 0.4*atom[t].sigma);
  }
  fclose(povfile);
  free(r);
}


int main(int argc, char *argv[])
{
  glas_system_t system;                 /* global system data */
  glas_atom_t *atom;                             /* atom data */
  int step, num_steps, 
    steps_eq;       /* number of steps used for equilibration */
  int argz=1;              /* counter for accessing arguments */
  int t, d;                                  /* loop counters */
  double time;                          /* time of simulation */
  double energy_pot, energy_kin;
  double T_start, T_end, factor;
  double virial0, virial1;    /* number, sum of virial values */
  double volume;
  int *n;
  

  system.dim = 3;                          /* initialisation */
  system.l = (double *) malloc(system.dim*sizeof(double));
  system.lh = (double *) malloc(system.dim*sizeof(double));
  n = (int *)  malloc(system.dim*sizeof(int));
  system.tau = 0.001;
  system.tau2 = system.tau*system.tau;
  system.T = 0.01;
  system.l[0] = 6.0*pow(2,1.0/6.0);
  system.l[1] = 0.5*system.l[0]*sqrt(3.0);
  system.l[2] = 0.5*system.l[0]*sqrt(8.0/3);
  volume = 1.0;
  for(d=0; d<system.dim; d++)
  {
    system.lh[d] = system.l[d]/2;
    volume *= system.l[d];
  }
  system.N = 6*6*6;
  n[0] = (int) floor(pow(system.N, 1.0/system.dim)+0.005);
  for(d=1; d<system.dim; d++)
    n[d] = n[0];
  sscanf(argv[argz++], "%lf", &system.T);
  num_steps = atoi(argv[argz++]);
  steps_eq = num_steps/2;
  factor = pow(T_end/T_start, 1.0/num_steps);
  atom = glas_setup(&system);
  glas_put_regular(&system, atom, n);
  
  virial0 = virial1 = 0;                 /* initialize virial */

  time = 0.0;                             /* start simulation */
  energy_pot = glas_energy_forces(&system, atom);
  energy_kin = glas_kinetic_energy(&system, atom);
  printf("%f %f %f %f %f %f\n", time, system.T, energy_pot/system.N, 
	 energy_kin/system.N, (energy_pot+energy_kin)/system.N,
	 system.virial);

  for(step=1; step<=num_steps; step++)           /* main loop */
  {
    energy_pot =  glas_velocity_verlet(&system, atom);
    glas_rescale_kinetic(&system, atom);
    energy_kin = glas_kinetic_energy(&system, atom);
    time += system.tau;
    if(step%100 == 0)
      printf("%f %f %f %f %f %f\n", time, system.T, 
	     energy_pot/system.N, energy_kin/system.N,
	     (energy_pot+energy_kin)/system.N, system.virial/volume);
    if(step>steps_eq)                       /* equilibrated ? */
    {
      virial1 += system.virial;
      virial0++;
    }
  }
  glas_plot_cfg(&system, atom);              /* make 3d image */

  printf("# T= %f, p= %f\n", system.T,
	 (system.N*system.T+virial1/virial0)/volume);

  free(system.l);                      /* clear up everything */
  free(system.lh);
  for(t=0; t<system.N; t++)
  {
    free(atom[t].x);
    free(atom[t].v);
    free(atom[t].f);
  }
  free(atom);
  free(n);
  return(0);
}