// Simple Molecular Dynamics Example Program

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

#include "Vector.h"
#include "Particles.h"
#include "Integrator.h"
#include "LennardJones.h"

// ------------------------------------------------------------

#define NR_PARTICLES 100 
#define NR_STEPS     10000
#define TIME_STEP    1.0
#define BOX_SIZE     12.0
#define MASS         5.0

// ------------------------------------------------------------

void SetupParticles( int N, Particle P[] )
{
 int i;

 srand( 0 );

 for( i = 0; i < N; ++i )
 {
  P[ i ].r.X = BOX_SIZE * ( rand() % 1000 - 500 ) / 500.0;
  P[ i ].r.Y = BOX_SIZE * ( rand() % 1000 - 500 ) / 500.0;
  P[ i ].r.Z = BOX_SIZE * ( rand() % 1000 - 500 ) / 500.0;

  P[ i ].p.X = rand() % 1000 - 500;
  P[ i ].p.Y = rand() % 1000 - 500;
  P[ i ].p.Z = rand() % 1000 - 500;

  P[ i ].p = ScalarVector( 5.0 / ( N * AbsVector( P[ i ].p ) ), P[ i ].p );

  P[ i ].m = MASS;
 }

 LennardJonesForces( N, P );
}

// ------------------------------------------------------------

int main( int argc, char *argv[] )
{
 Particle *P1, *P2, *PP;
 int i;

 printf( "Simple Molecular Dynamics Example Program\n\n" );

 P1 = calloc( NR_PARTICLES, sizeof( Particle ) );
 P2 = calloc( NR_PARTICLES, sizeof( Particle ) );

 printf( "Setup of the particles...\n\n" );

 SetupParticles( NR_PARTICLES, P1 );

 for( i = 0; i < NR_PARTICLES; ++i )
  PrintParticle( P1[ i ] );

 printf( "\n\nStart of dynamics run...\n\n" );

 for( i = 0; i < NR_STEPS; ++i )
 {
  IntegrateATimeStep( NR_PARTICLES, P1, P2, TIME_STEP, LennardJonesForces );

  if( !( i % ( NR_STEPS / 100 )  ) ) printf( "t: %lE, Ek: %lE, V: %lE, E: %lE\n", ( TIME_STEP * i ), TotalKineticEnergy( NR_PARTICLES, P1 ),
                                        LennardJonesPotential( NR_PARTICLES, P1 ), TotalKineticEnergy( NR_PARTICLES, P1 ) +
                                        LennardJonesPotential( NR_PARTICLES, P1 ) );

  PP = P2;
  P2 = P1;
  P1 = PP;
 }

 free( P1 );
 free( P2 );
 
 printf( "\n\nDone...\n\n" ); 
}