/*
  simctop.c

  C version of a simulation of a fluid of molecules each consisting of
  four Lennard-Jones sites.

  This program is meant only as an illustration of the use of the
  implementation of the exact rotation of a rigid body in the context
  of the symplectic integration scheme using the exact rotation matrix
  (SIER2). (The calculations of the forces and the random number
  generator is therefore implemented only in the most straightforward
  and inefficient way.)

  Compile the program with
    cc simtop.c ctops.c -lgsl -lgslcblas -lm -o simtop

  Note that ctops.h requires the GNU Scientific Libray to be installed
  and linked.

  Run with "simtop <inifile>" or just "simtop" to use "simtop.ini" as
  an inifile (initialization file). 

  This inifile should contain a list of numbers which signify

     number of particles
     linear size of the cubic simulation box
     run-time of the simulation
     integration step-size 

  The program writes a line per time step containing four columns representing
  time, total energy, translational kinetic energy, rotational kinetic
  energy and potential energy. It is written to the standard output.

  Ramses van Zon, May 11, 2007
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "ctops.h" /* functions computing exact rotation in free motion */

/* macros for some Vector and Matrices manipulations */
#define VCROSS(c,a,b) ((c).x=(a).y*(b).z-(a).z*(b).y,\
                       (c).y=(a).z*(b).x-(a).x*(b).z,\
                       (c).z=(a).x*(b).y-(a).y*(b).x)
#define MVMUL(c,M,a)  ((c).x=(M).xx*(a).x+(M).xy*(a).y+(M).xz*(a).z,\
                       (c).y=(M).yx*(a).x+(M).yy*(a).y+(M).yz*(a).z,\
                       (c).z=(M).zx*(a).x+(M).zy*(a).y+(M).zz*(a).z)
#define MTVMUL(c,M,a) ((c).x=(M).xx*(a).x+(M).yx*(a).y+(M).zx*(a).z,\
                       (c).y=(M).xy*(a).x+(M).yy*(a).y+(M).zy*(a).z,\
                       (c).z=(M).xz*(a).x+(M).yz*(a).y+(M).zz*(a).z)

/* random number generator*/
double rnd(void)
{
  return rand()/(double)RAND_MAX; 
}

/* molecular properties */
struct Molecule {
  Vector q;  /* position */
  Vector p;  /* momentum */
  Matrix A;  /* attitude */
  Vector L;  /* angular momentum */
  Vector F;  /* force */
  Vector t;  /* torque */
  double V;  /* potential energy */
};
typedef struct Molecule Molecule;

Top top;     /* see ctops.h and docctops.pdf */

#define nSites 4
Vector site[nSites] = {{     0,  0.075, 0.05},    /* sites in the molecule */
	 	       {     0, -0.075, 0.05},
		       { 0.075,      0, -0.1},
		       {-0.075,      0, -0.1} };
Vector I = {0.3375,0.45,0.525}; /* moments of inertia of a molecule */
double m = 60.0;                /* mass of a molecule */

int N;              /* number of particles */
double L;           /* linear length of the (cubic) system */
double RUNTIME;     /* runtime */
double DT;          /* integration time step */

Molecule *mol;      /* the molecules */

/* 
   The site-site potential is cut-off as follows:
     for r < RCA, full Lennard-Jones is used;
     for r > RCB, the potential and force are both zero;
     for RCA < r < RCB, an interpolating cubic is used.
   As a result, the potential nor the forces have any discontinuities.
*/
#define RCA 2.5
#define RCB 4
#define RCA2 6.25
#define RCB2 16
double cubic[4];  /* coefficients in the cubic interpolation */
void computeCubic(void)
{
  double ir2, ir6, V, F, a, b, c, ia2, ia3;

  ir2 = 1/(double)RCA2;
  ir6 = ir2 * ir2 * ir2;
  V   = ir6 * (4 * ir6 - 4);
  F   = ir6 * (48 * ir6 - 24) / RCA;
  a   = RCB - RCA;
  b   = 3 * V - a * F;
  c   = -(2 * V - a * F);
  ia2 = 1/(a*a);
  ia3 = ia2/a;
  cubic[0]= b * RCB2 * ia2 + c * RCB * RCB2 * ia3;
  cubic[1]= -2 * b * RCB * ia2 - 3 * c * RCB2 * ia3;
  cubic[2]= b * ia2 + 3 * c * RCB * ia3;
  cubic[3]= -c * ia3;
}

void siteForce(Vector *r, Vector *f, double *V)
{
  double rnrm, r2, ir2, ir6, F;

  r2 = r->x * r->x +  r->y * r->y + r->z * r->z;

  if (r2 > RCB2) { 
     *V = 0.0; 
     f->x = f->y = f->z = 0;
  } else if (r2 > RCA2) { 
    /* smooth interpolation */
    rnrm  = sqrt(r2);
    *V = cubic[0] + rnrm * cubic[1] + r2 * (cubic[2] + rnrm * cubic[3]);
    F  = -(cubic[1] / rnrm + 2 * cubic[2] + 3 * rnrm * cubic[3]);
    f->x = F * r->x;
    f->y = F * r->y;
    f->z = F * r->z;
  } else { 
    ir2 = 1/r2;
    ir6 = ir2 * ir2 * ir2;
    *V  = ir6 * (4 * ir6 - 4);
    F  = ir6 * ir2 * (48 * ir6 - 24);
    f->x = F * r->x;
    f->y = F * r->y;
    f->z = F * r->z;
  } 
}

/* force F and torque t on molecule B due to molecule A.
   r should be the distance vector from the center of mass of A 
   to that of B. */
void molForce(Molecule *A, Molecule *B, Vector *r,
	      Vector *F, Vector *t, double *V)
{
  int a, b;   
  double vba; 
  Vector rba, fba, tba, ra, rb;
  
  *V = 0;
  F->x = F->y = F->z = 0;
  t->x = t->y = t->z = 0;

  for (a = 0; a < nSites; a++) {
    for (b = 0; b < nSites; b++) { 
      /* construct vector pointing from a to b */ 
      MTVMUL ( ra, A->A, site[a] );
      MTVMUL ( rb, B->A, site[b] );

      rba.x = r->x + rb.x - ra.x;
      rba.y = r->y + rb.y - ra.y;
      rba.z = r->z + rb.z - ra.z;

      siteForce ( &rba, &fba, &vba );

      *V += vba;
      F->x += fba.x;
      F->y += fba.y;
      F->z += fba.z;
      VCROSS ( tba, rb, fba );

      t->x += tba.x; 
      t->y += tba.y;
      t->z += tba.z;  
    }
  }
}

void computeForces(void) /* straightforward but slow */
{
  int i, j;    
  double V;    
  Vector r;    
  Vector F;    
  Vector t;    
  Vector tcm;  /* torque due to different centers of mass */
  
  for (i = 0; i < N; i++){
    mol[i].F.x = mol[i].F.y = mol[i].F.z = 0;
    mol[i].t.x = mol[i].t.y = mol[i].t.z = 0;
    mol[i].V = 0;
  }

  for (i = 0; i < N; i++) {
    for (j = 0; j < i; j++) {
      /* determine distance of j from i */
      r.x = mol[j].q.x - mol[i].q.x;  
      r.y = mol[j].q.y - mol[i].q.y;  
      r.z = mol[j].q.z - mol[i].q.z;  

      /* apply periodic boundary conditions */
      r.x -= L*rint(r.x/L);
      r.y -= L*rint(r.y/L);
      r.z -= L*rint(r.z/L);
      
      /* compute force and torque on j due to i */
      molForce( &(mol[i]), &(mol[j]), &r, &F, &t, &V );
      
      mol[j].F.x += F.x;
      mol[j].F.y += F.y;
      mol[j].F.z += F.z;
      mol[i].F.x -= F.x;
      mol[i].F.y -= F.y;
      mol[i].F.z -= F.z;
      mol[j].t.x += t.x;
      mol[j].t.y += t.y;
      mol[j].t.z += t.z;
      /* the torque on i is minus the torque on j, minus r x Fi: */
      VCROSS   ( tcm, r, F );
      mol[i].t.x -= t.x + tcm.x;
      mol[i].t.y -= t.y + tcm.y;
      mol[i].t.z -= t.z + tcm.z;

      mol[i].V += 0.5 * V;
      mol[j].V += 0.5 * V;
    }
  }
}

/* read in parameters from a file and initialize the molecules */
void initialize(char *filename)
{
  FILE *f; 
  int i;

  f = fopen(filename, "r");
  if (f == 0) {
    fprintf(stderr, "Cannot find file %s\nExiting.\n", filename);
    exit(1);
  }

  fscanf(f, "%d", &N);
  fscanf(f, "%lf", &L);
  fscanf(f, "%lf", &RUNTIME);
  fscanf(f, "%lf", &DT);

  DefineTop( &top, RecurTop, I.x, I.y, I.z );

  mol = malloc(N * sizeof (Molecule));

  for (i = 0; i < N; i++) {
    /* random initial positions */
    mol[i].q.x = L * (rnd() - 0.5); /* positions range from -L/2 to L/2 */
    mol[i].q.y = L * (rnd() - 0.5);
    mol[i].q.z = L * (rnd() - 0.5);

    /* molecule initially in upright attitude */
    mol[i].A.xx = 1.0;  mol[i].A.xy = 0.0;  mol[i].A.xz = 0.0;
    mol[i].A.yx = 0.0;  mol[i].A.yy = 1.0;  mol[i].A.yz = 0.0;
    mol[i].A.zx = 0.0;  mol[i].A.zy = 0.0;  mol[i].A.zz = 1.0;

    /* random initial momenta */
    mol[i].p.x = sqrt(m) * (rnd() - 0.5);
    mol[i].p.y = sqrt(m) * (rnd() - 0.5);
    mol[i].p.z = sqrt(m) * (rnd() - 0.5);

    /* random initial angular momenta */
    mol[i].L.x = sqrt(I.x) * (rnd() - 0.5);
    mol[i].L.y = sqrt(I.y) * (rnd() - 0.5);
    mol[i].L.z = sqrt(I.z) * (rnd() - 0.5);
  }
  computeCubic();
  computeForces();
}

/* take a single time-step */
void timeStep(void)
{
  int i;
  Vector omega;

  /* force propagatation by half a step */
  for (i = 0; i < N; i++) {
    mol[i].p.x += 0.5 * DT * mol[i].F.x;
    mol[i].p.y += 0.5 * DT * mol[i].F.y;
    mol[i].p.z += 0.5 * DT * mol[i].F.z;
    mol[i].L.x += 0.5 * DT * mol[i].t.x;
    mol[i].L.y += 0.5 * DT * mol[i].t.y;
    mol[i].L.z += 0.5 * DT * mol[i].t.z;
  }

  /* free propagation by a full step */
  for ( i = 0; i < N; i++ ) {
    MVMUL ( omega, mol[i].A, mol[i].L ); 
    omega.x /= I.x;
    omega.y /= I.y; 
    omega.z /= I.z; 
    Propagation( &top, DT, &omega, &(mol[i].A) ); 

    mol[i].q.x += DT / m * mol[i].p.x;
    mol[i].q.y += DT / m * mol[i].p.y;
    mol[i].q.z += DT / m * mol[i].p.z;
  }

  /* recompute forces */
  computeForces();

  /* force propagatation by half a step */
  for (i = 0; i < N; i++) {
    mol[i].p.x += 0.5 * DT * mol[i].F.x;
    mol[i].p.y += 0.5 * DT * mol[i].F.y;
    mol[i].p.z += 0.5 * DT * mol[i].F.z;
    mol[i].L.x += 0.5 * DT * mol[i].t.x;
    mol[i].L.y += 0.5 * DT * mol[i].t.y;
    mol[i].L.z += 0.5 * DT * mol[i].t.z;
  }
}

/* report energetic properties to standard output */
void report(double t)
{
  int i;
  double Etrans, Erot, Epot;
  Vector Lb;
  
  Etrans = Erot = Epot = 0;
  for (i = 0; i < N; i++) {
    MVMUL ( Lb, mol[i].A, mol[i].L );
    Epot += mol[i].V;
    Etrans += 0.5 * (mol[i].p.x * mol[i].p.x 
		     + mol[i].p.y * mol[i].p.y 
		     + mol[i].p.z * mol[i].p.z) / m;
    Erot += 0.5*( Lb.x * Lb.x / I.x 
                + Lb.y * Lb.y / I.y 
                + Lb.z * Lb.z / I.z);
  }
  printf("%16.14e %16.14e %16.14e %16.14e %16.14e\n", 
	 t, Etrans+Erot+Epot, Etrans, Erot, Epot);
}

/* run the simulation initialized by initialize(...) */
void run(void)
{
  double t;
  
  for (t = 0; t < RUNTIME; t += DT) {
    report(t);
    timeStep();
  }
  report(t);
}


int main(int argc, char **argv)
{
  char def_ini[] = "simtop.ini";

  if (argc > 1)    
    initialize(argv[1]);
  else 
    initialize(def_ini);

  run();

  return 0;
}