//
// simtop.cc
//
// 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 are therefore implemented only in the most
// straightforward and inefficient way.)
//
// Compile the program with
//   c++ simtop.cc tops.cc -lgsl -lgslcblas -lm -o simtop
//
// Note that tops.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 five 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 <iostream>
#include <iomanip>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include "tops.h" // functions computing exact rotation in free motion 
using namespace std; 

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

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

class System
{
  static const int nSites = 4;
  static const Vector site[nSites];  // sites in the molecule 
  static const Vector I;    // moments of inertia of a molecule
  static const double m;    // mass of a molecule
  static double cubic[4];   // coefficients in the cubic interpolation 
  static TopRecur top;      // see tops.h and doctops.pdf

  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
  void computeCubicCoefficients()
  {
    double ir2 = 1/(double)RCA2;
    double ir6 = ir2 * ir2 * ir2;
    double V   = ir6 * (4 * ir6 - 4);
    double F   = ir6 * (48 * ir6 - 24) / RCA;
    double a   = RCB - RCA;
    double b   = 3 * V - a * F;
    double c   = -(2 * V - a * F);
    double ia2 = 1/(a*a);
    double 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(const Vector &r, Vector &f, double &V) const
  {
    double r2 = r.nrm2();

    if (r2 > RCB2) { 
      V = 0.0; 
      f.zero(); 
    } else if (r2 > RCA2) { 
      // smooth interpolation 
      double 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]) * r;
    } else { 
      double ir2 = 1/r2;
      double ir6 = ir2 * ir2 * ir2;
      V  = ir6 * (4 * ir6 - 4);
      f  = ir6 * ir2 * (48 * ir6 - 24) * r;
    } 
  }
  
  // Compute the force F and torque t on molecule B due to molecule A.
  // Note: r should be the distance vector from the center of mass of
  // A to that of B.
  void molForce(const Molecule &A, const Molecule &B, const Vector &r,
		Vector &F, Vector &t, double &V) const
  {
    V = 0;          
    F.zero();
    t.zero(); 

    for (int a = 0; a < nSites; a++) {
      for (int b = 0; b < nSites; b++) { 
	Vector ra = Transpose(A.A) * site[a];
	Vector rb = Transpose(B.A) * site[b];
	Vector rba = r + rb - ra;
	Vector fba;
	double vba;
	
	siteForce ( rba, fba, vba );
	
	V += vba;
	F += fba;
	t += rb^fba;
      }
    }
  }

  void computeForces()   // straightforward but slow 
  {
    for (int i = 0; i < N; i++){
      mol[i].F.zero();
      mol[i].t.zero();
      mol[i].V = 0;
    }
    
    for (int i = 0; i < N; i++) {
      for (int j = 0; j < i; j++) {
	Vector r =  mol[j].q - mol[i].q;  
	
	// 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 
	double V;    
	Vector F;
	Vector t;
	molForce( mol[i], mol[j], r, F, t, V );
	
	mol[j].F += F;
	mol[i].F -= F;
	mol[j].t += t;
	mol[i].t -= t + (r^F); // torque on i is -torque on j, - r x Fj: 
      
	mol[i].V += 0.5 * V;
	mol[j].V += 0.5 * V;
      }
    }
  }
  
  // take a single time-step 
  void timeStep()
  {
    // force propagatation by half a step 
    for (int i = 0; i < N; i++) {
      mol[i].p += 0.5 * DT * mol[i].F;
      mol[i].L += 0.5 * DT * mol[i].t;
    }

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

    // recompute forces 
    computeForces();

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

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

 public:

  System(): mol(0)
  {
    cout << setprecision(14);   // cause accurate output
    computeCubicCoefficients(); // precompute coefficients
  }

  ~System()
  {
    delete[] mol;
  }

  void initialize(const char *filename)
  {
    // read in parameters 
    ifstream f(filename); 
    if (!f) {
      cerr << "Cannot find file " << filename << "\nExiting.\n";
      exit(1);
    }
    f >> N;
    f >> L;
    f >> RUNTIME;
    f >> DT;

    //  initialize the system 
    mol = new Molecule[N];

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

      mol[i].A.one();         // molecule initially in upright attitude 

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

      mol[i].L.x = sqrt(I.x)*(rnd()-0.5); // random angular momenta 
      mol[i].L.y = sqrt(I.y)*(rnd()-0.5); // scaled with respective
      mol[i].L.z = sqrt(I.z)*(rnd()-0.5); // moments of inertia
    }

    computeForces();
  }

  void run()
  {
    report(0);
    for (double t = 0; t < RUNTIME; t += DT) {
      timeStep();
      report(t);
    }
  }
};

const Vector System::site[nSites] 
   = {Vector(    0,  .075,  .05),
      Vector(    0, -.075,  .05),
      Vector( .075,     0, -.1),
      Vector(-.075,     0, -.1) };
const Vector System::I(0.3375,0.45,0.525);
const double System::m = 60.0;            
TopRecur System::top(System::I.x, System::I.y, System::I.z);
double System::cubic[4];


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

  System system;

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

  system.run();

  return 0;
}