//code for tool used to calculate the initial potential energy of a block of particles used initialized by the Verlet molecular dynamics code.
#include <vector.h>
#include <stdio.h>
#include <math.h>
#include <fstream.h>

class Particle {
   public:
      Particle();
      Particle( const double x, const double y, const double z ) {
	m_x = x; 
	m_y = y;
	m_z = z;
	m_potential = 0.0;
      }
      ~Particle() {
      };
      double Getx() {
	 return m_x;
      }
      double Gety() {
	 return m_y;
      }
      double Getz() {
	 return m_z;
      }
      void updatePotential ( const double potential ) {
	 m_potential += potential;
      }
      double GetPotential() {
	 return m_potential;
      }
   private:
      double m_x;
      double m_y;
      double m_z;
      double m_potential;
};

int main () {

   int numberOfParticles = 0;
   double density = 0.0;
   printf ( "Enter number of particles.\n N = " );
   scanf( "%d", &numberOfParticles );
   printf( "Enter the density.\n rho = " );
   scanf( "%lf", &density );

   double bigEdge = pow ( numberOfParticles/density, 1./3. );
   double littleEdge = bigEdge/pow( ((double)(numberOfParticles)), 1./3. );
   
   vector<Particle> particleVector;

   for ( int unsigned i=0; i<((int unsigned)(rint(bigEdge/littleEdge))); i++ ) {
      for ( int unsigned j=0; j<((int unsigned)(rint(bigEdge/littleEdge))); j++ ) {
	for ( int unsigned k=0; k<((int unsigned)(rint(bigEdge/littleEdge))); k++ ) {
	  Particle aParticle( -bigEdge/2.+littleEdge/2.+((double)(k))*littleEdge, -bigEdge/2.+littleEdge/2.+((double)(j))*littleEdge, -bigEdge/2.+littleEdge/2.+((double)(i))*littleEdge );
	  particleVector.push_back( aParticle );  
	 }
      }
   }

   printf( "%d", particleVector.size() );
   for ( vector<Particle>::iterator particle = particleVector.begin(); particle != particleVector.end(); particle++ ) {
      printf( "x = %lf, y = %lf, z = %lf\n", particle->Getx(), particle->Gety(), particle->Getz() );
   }
   
   double separationSquared = 0.0;
   double separationx = 0.0;
   double separationy = 0.0;
   double separationz = 0.0;

   ofstream tempFile( "Tmp.dat", ios::app );

   for ( vector<Particle>::iterator i = particleVector.begin(); i != particleVector.end()-1; i++ ) {
      for ( vector<Particle>::iterator j = i+1; j != particleVector.end(); j++ ) {
	separationx = (i->Getx()-j->Getx());
	separationx -= bigEdge*rint(separationx/bigEdge);
	separationy = (i->Gety()-j->Gety());
	separationy -= bigEdge*rint(separationy/bigEdge);
	separationz = (i->Getz()-j->Getz());
	separationz -= bigEdge*rint(separationz/bigEdge);

	separationSquared = (separationx*separationx)+(separationy*separationy)+(separationz*separationz);
	//tempFile << separationSquared << endl;
	i->updatePotential( 4*(pow( separationSquared, -6.0 )-pow( separationSquared, -3.0 ) ) );
	j->updatePotential( 4*(pow( separationSquared, -6.0 )-pow( separationSquared, -3.0 ) ) );
	tempFile << i->GetPotential() << "  " << j->GetPotential() << endl;
     }
   }

   double averagePotential = 0.0;

   for ( vector<Particle>::iterator particle = particleVector.begin(); particle != particleVector.end(); particle++ ) {
       averagePotential += particle->GetPotential();
   }
   averagePotential /= ((double)(numberOfParticles));
   averagePotential -= 4*(pow( bigEdge/2., -12.0 )-pow( bigEdge/2., -6.0 ));

   printf ( "Initial Potential Per Particle = %lf \n", averagePotential );

   tempFile.close();

}