/*   []----------------------------------------[]
      |                Main.cpp                |
     []----------------------------------------[]
      |                                        |
      |   AUTHOR:    MFSomers 2005.            |
      |   USE:       Main program ...          |
      |                                        |
     []----------------------------------------[]
*/

// Copyright (C) 2005 M.F. Somers, Theoretical Chemistry Department, Leiden University
//
// This is free software; you can redistribute it and/or modify it under the terms of
// the GNU General Public License as published by the Free Software Foundation.
//
// http://www.gnu.org/licenses/gpl.txt


// Specify your compile options here if you want to... you can also put them in the makefile...

// #define INPUT_FILE_CHECKER
// #define USE_BOINC_API
// #define USE_GLUT_GRAPH_API
// #define USE_BOINC_GRAPH_API

#if defined(INPUT_FILE_CHECKER)

#   if defined(USE_BOINC_API) || defined(USE_GLUT_GRAPH_API) || defined(USE_BOINC_GRAPH_API)
#      error Cannot use any of the API for the input file check program ...
#   endif

#else

#if defined(USE_BOINC_API) && defined(USE_GLUT_GRAPH_API)
#   error Cannot use GLUT API graphics with the non-graphical BOINC API ...
#endif

#if defined(USE_GLUT_GRAPH_API) && defined(USE_BOINC_GRAPH_API)
#   error Either user the GLUT API or the BOINC graphical API ...
#endif

// use boinc automatically when using graphical boinc
#if defined(USE_BOINC_GRAPH_API)
#   define USE_BOINC_API
#endif

#endif

// first include system specific things 
#if defined(__linux__)
#   if defined(USE_BOINC_GRAPH_API)
#      define _REENTRANT               // aparently this is needed according to yello example
#   endif
#elif defined(_WIN32) || defined(WIN32)
#elif defined(__APPLE__)
#elif defined(FREEBSD)
#else
#     error This platform is not suported !
#endif

// default includes we need on all systems or platforms
#include <ctime>

#if defined(USE_CRLIBM)
#include "crlibm.h"
#endif

#if defined(USE_BOINC_API)          // so we can use them here
void CreateCheckPointFile( void );
#endif

/* ------------------------------------------------ My general includes ---- */

#include "Global.h"
#include "Vector.h"
#include "Rand.h"
#include "Constants.h"
#include "ClassicalDynamics.h"

#include "PointPair.h"
#include "Coulomb.h"
#include "Gravity.h"
#include "Harmonic.h"
#include "Morse.h"
#include "LennardJones.h"
#include "Rydberg.h"

#include "NewtonDynamics.h"
#include "LangevinDynamics.h"
#include "SteepestDescent.h"

#include "Bending.h"
#include "Torsional.h"

#include "ReadInput.h"

/* --------------------------------------------- My global defines ------- */

#define ENERGY_CONSERVATION_CHECK      1.0E-06
#define MAX_TIME_DELTA                 0.01
#define END_TIME                       100.0
#define BOX_SIZE                       25.0
#define MAX_CONFORMATION_STEP_SIZE     1.0
#define MAX_CONFORMATION_NSTEPS        10000
#define CONFORMATION_DEFAULT_ACCURACY  1.0E-06

const Vector ONE_VECTOR = { 1.0, 1.0, 1.0 };

/* -------------------------------------- My newton dynamics class ------- */

class MyLangevinDynamics;
typedef MyLangevinDynamics *MyLangevinDynamicsPointer;

class MyLangevinDynamics : public LangevinDynamics
                         {
			  public:
			  
			  MyLangevinDynamics( double TheStartTime = 0.0, double maxdt = MAX_TIME_DELTA, double tend = END_TIME, 
                                              Vector TheBoxVector = ONE_VECTOR, double dE = ENERGY_CONSERVATION_CHECK, 
                                              int ThePeriodicity = 0, int UseConstantTemperature = 0, double TheT = 0.0,
					      double Bolzmann = 0.0, double TheTReScaleTime = 0.0, double RandomMFactor = 0.0, 
                                              double TheGamma = 0.0, double TheCMRMF = 0.0, unsigned int TheId = GET_NEW_ID ) : LangevinDynamics( TheGamma, Bolzmann, TheT, TheId )
			  {
			   Time = StartTime = TheStartTime;                           
			   NewTimeStep = MaxTimeStep = TimeStep = ( maxdt < CONST_EPSILON ? CONST_EPSILON : maxdt );
			   EndTime = ( tend > StartTime + TimeStep ? tend : StartTime + TimeStep );
			   PreviousTime = StartTime - TimeStep;
			   
			   if( TheBoxVector.X < CONST_EPSILON ) TheBoxVector.X = CONST_EPSILON;
			   if( TheBoxVector.Y < CONST_EPSILON ) TheBoxVector.Y = CONST_EPSILON;
			   if( TheBoxVector.Z < CONST_EPSILON ) TheBoxVector.Z = CONST_EPSILON;
			   TheBox = TheBoxVector;
			   
			   Periodicity = ThePeriodicity;
			   
			   // ConstantTemperature -1 means initialize at t=0, 0 means do nothing, >= 1 means rescale velocities and stuff like that
			   ConstantTemperature = UseConstantTemperature;
			   TargetTemperature = ( TheT < CONST_EPSILON ? ConstantTemperature = 0, 0.0 : TheT );
			   BoltzmannConstant = ( Bolzmann < CONST_EPSILON ? ConstantTemperature = 0, 0.0 : Bolzmann );
			   TemperatureReScaleTime = ( TemperatureReScaleTime < CONST_EPSILON ? 0.0 : TemperatureReScaleTime );
			   if( ConstantTemperature & 2 ) TemperatureReScaleTime = ( TheTReScaleTime >= MaxTimeStep ? TheTReScaleTime : MaxTimeStep );
			   RandomMomentaFactor = ( RandomMFactor < CONST_EPSILON ? 0.0 : RandomMFactor );
                           CenterOfMassRemovalFactor = ( TheCMRMF < CONST_EPSILON ? 0.0 : TheCMRMF );
			   			   
			   Accuracy = ( dE < CONST_EPSILON ? CONST_EPSILON : dE );
			   			   
			   TotalNrOfStepsDone = 0;
			   TheTotalKineticEnergy = TheTotalInteractionPotential = 0.0;
			   TheAccumulatedTotalKineticEnergy = TheAccumulatedTotalInteractionPotential = 0.0;
			   TheTemperature = ThePressure = 0.0;
			   TheAccumulatedTemperature = TheAccumulatedPressure = 0.0;
			   
			   TheCentreOfMass.X = TheCentreOfMass.Y = TheCentreOfMass.Z = 0.0;
			   TheDipoleMoment.X = TheDipoleMoment.Y = TheDipoleMoment.Z = 0.0;
			   		   
			   NrOfStepsToDo = 1; NrOfTimeStepsDone = 0; FindAppropriateTimeStep = 1; StopTrajectory = 0;
			   
			   if( ConstantTemperature > 0 ) FindAppropriateTimeStep = 0;
			  }
			  
			  virtual int IntegratorConservationCheck( double Step );
			  virtual int TrajectoryContinuationCheck( double At, double Step, int nStep, int nSteps );
			  
			  int          Periodicity, ConstantTemperature;
                          Vector       TheCentreOfMass, TheDipoleMoment, TheBox;
			  
			  double       PreviousTime, StartTime, Time, TimeStep, NewTimeStep, MaxTimeStep, EndTime;
                          int          NrOfStepsToDo, NrOfTimeStepsDone, FindAppropriateTimeStep, StopTrajectory;
			  
			  double       Accuracy;
			  		  
			  unsigned int TotalNrOfStepsDone;
			  double       TheTotalKineticEnergy, TheTotalInteractionPotential, RandomMomentaFactor, CenterOfMassRemovalFactor;
                          double       TheAccumulatedTotalKineticEnergy, TheAccumulatedTotalInteractionPotential;
			  double       TargetTemperature, BoltzmannConstant, TheTemperature, ThePressure;
			  double       TemperatureReScaleTime, TheAccumulatedTemperature, TheAccumulatedPressure;
			 };

/* -------------------------------------- Global data of this program ---- */

char                    OutputFileName[ 256 ],
                        InputFileName[ 256 ],
                        StdOutFileName[ 256 ],
                        StdErrFileName[ 256 ];
FILE                    *InputFilePointer = NULL,
                        *OutputFilePointer = NULL,
                        *StdOutFilePointer = NULL,
                        *StdErrFilePointer = NULL;
SteepestDescentPointer  Conformation = NULL;
MyLangevinDynamicsPointer Dynamics = NULL;
InputFileDataType       TheInputData;

/* --------------------------------------- The newtondynamics methods ---- */
/* The conservation of energy integration check                            */

int MyLangevinDynamics::IntegratorConservationCheck( double Step )
{
 double OldEnergy, NewEnergy;
 
 if( FindAppropriateTimeStep && ( ConstantTemperature <= 0 ) )
 {
     OldEnergy = TotalPotentialEnergy( *GetStateInteractionsPointer( 1 ) ) + 
                 TotalKineticEnergy( *GetStateParticlesPointer( 1 ) );

     NewEnergy = TotalPotentialEnergy( *GetStateInteractionsPointer( 0 ) ) + 
                 TotalKineticEnergy( *GetStateParticlesPointer( 0 ) );	     

     if( fabs( OldEnergy - NewEnergy ) >= Accuracy * fabs( OldEnergy ) )    // then check if energy is conserved
      return( 0 ); 
    
     NewTimeStep = Step;                               // if so use this step size in future runs
     
     FindAppropriateTimeStep = 0;
 }
  
 return( 1 );
}

/* ----------------------------------------------------------------------- */
/* The trajectory continuation check to get optimal timesteps and so on    */

int MyLangevinDynamics::TrajectoryContinuationCheck( double At, double Step, int nStep, int nSteps )
{
 ParticleListPointer TheState = GetFinalStateParticlesPointer();
 ParticlePointer     TheParticle;
 Vector              ThePosition, CenterOfMassVelocity;
 double              MaxXValue, MaxYValue, MaxZValue, ScaleFactor, TotalMass;
 int                 theNrOfParticlesIncluded = NrOfParticlesIncluded();
 int                 Size = TheState -> size(); 
 
 Time = At;

 if( ( nStep >= nSteps ) || ( NewTimeStep < Step ) )      // check if step has changed or should not continue with the run
 {
  // calculate some properties at the end of each trajectory run

  TheCentreOfMass              = CentreOfMass( (*TheState) );  
  TheDipoleMoment              = DipoleMoment( (*TheState) );      
  TheTotalKineticEnergy        = TotalKineticEnergy( (*TheState) );
  TheTotalInteractionPotential = TotalPotentialEnergy( (*GetFinalStateInteractionsPointer()) ); 
  
  TheTemperature = ( BoltzmannConstant > 0.0 ? 0.66666666666666 * TheTotalKineticEnergy / ( theNrOfParticlesIncluded * BoltzmannConstant ) : 0.0 );
  ThePressure = ( 0.66666666666666 * TheTotalKineticEnergy + 0.3333333333333 * Virial( (*TheState), TheBox ) ) / ( TheBox.X * TheBox.Y * TheBox.Z );
  
  TheAccumulatedTotalKineticEnergy += TheTotalKineticEnergy;
  TheAccumulatedTotalInteractionPotential += TheTotalInteractionPotential;  
  TheAccumulatedTemperature += TheTemperature;
  TheAccumulatedPressure += ThePressure;
  ++TotalNrOfStepsDone;
  
  if( !Periodicity )
  {
   
   // and check if any of the particle runs out of the box, if we do not have periodicity
   for( int n = 0; n < Size; ++n )
    if( (*(TheState))[ n ] != NULL )
    {
     ThePosition = (*(TheState))[ n ] -> Position;
     
     MaxXValue = 0.5 * TheBox.X + 1.25 * (*(TheState))[ n ] -> Radius;
     MaxYValue = 0.5 * TheBox.Y + 1.25 * (*(TheState))[ n ] -> Radius;
     MaxZValue = 0.5 * TheBox.Z + 1.25 * (*(TheState))[ n ] -> Radius;
     
     if( ThePosition.X >= MaxXValue ) StopTrajectory = 1;
     if( ThePosition.X <= -MaxXValue ) StopTrajectory = 1;
     if( ThePosition.Y >= MaxYValue ) StopTrajectory = 1;
     if( ThePosition.Y <= -MaxYValue ) StopTrajectory = 1;    
     if( ThePosition.Z >= MaxZValue ) StopTrajectory = 1;
     if( ThePosition.Z <= -MaxZValue ) StopTrajectory = 1;    
    }
    
  }
  
  // we should rescale the velocities to a target temperature or stuff like that...
  if( ConstantTemperature > 0 )                  
  {
  
   // check if we want to add some randomness to the momenta... only if not doing langevin dynamics...
   if( ( RandomMomentaFactor >= CONST_EPSILON  ) && ( Gamma < CONST_EPSILON ) && ( ConstantTemperature & 4 ) && !( ConstantTemperature & 16 ) )
   {
    // re-initialize some data...
    TheTotalKineticEnergy = 0.0;
    TheAccumulatedTemperature -= TheTemperature;

    // add some random momenta to all particles...
    for( int n = 0; n < Size; ++n )
     if( ( TheParticle = (*(TheState))[ n ] ) != NULL )
      if( TheParticle -> Mass > 0.0 )
      {
       TheParticle -> Momentum += RandomMomentaFactor * Abs( TheParticle -> Momentum ) * ARandomNormalizedVector();
       TheTotalKineticEnergy += 0.5 * Dot( TheParticle -> Momentum, TheParticle -> Momentum ) / TheParticle -> Mass;
      }

    // finaly re-calculate the temperature... 
    TheTemperature = ( BoltzmannConstant > 0.0 ? 0.66666666666666 * TheTotalKineticEnergy / ( theNrOfParticlesIncluded * BoltzmannConstant ) : 0.0 );
    TheAccumulatedTemperature += TheTemperature;
   }



   // check if we want to remove center of mass velocity...
   if( ( CenterOfMassRemovalFactor >= CONST_EPSILON ) && ( ConstantTemperature & 8 ) )
   {
    // re-initialize some data...
    TheAccumulatedTemperature -= TheTemperature;
    TheTotalKineticEnergy = 0.0;
    TotalMass = 0.0;
    CenterOfMassVelocity.X = CenterOfMassVelocity.Y = CenterOfMassVelocity.Z = 0.0;

    // now calculate total center of mass velocity
    for( int n = 0; n < Size; ++n )
     if( ( TheParticle = (*(TheState))[ n ] ) != NULL )
     {
      TotalMass += TheParticle -> Mass;
      CenterOfMassVelocity += TheParticle -> Momentum;
     }

    if( TotalMass > 0.0 )
     CenterOfMassVelocity /= TotalMass;
    else
     CenterOfMassVelocity *= 0.0;

    // take the set factor into account...
    CenterOfMassVelocity *= CenterOfMassRemovalFactor;

    // remove total center of mass velocity from all particles to get total momentum
    // to be zero and then calculate the current total kinetic energy of the system
    for( int n = 0; n < Size; ++n )
     if( ( TheParticle = (*(TheState))[ n ] ) != NULL )
      if( TheParticle -> Mass > 0.0 )
      {
       TheParticle -> Momentum -= CenterOfMassVelocity * TheParticle -> Mass;
       TheTotalKineticEnergy += 0.5 * Dot( TheParticle -> Momentum, TheParticle -> Momentum ) / TheParticle -> Mass;
      }
 
    // finaly re-calculate the temperature... 
    TheTemperature = ( BoltzmannConstant > 0.0 ? 0.66666666666666 * TheTotalKineticEnergy / ( theNrOfParticlesIncluded * BoltzmannConstant ) : 0.0 );
    TheAccumulatedTemperature += TheTemperature;
   }




   // do rescaling if needed unless when doing langevin dynamics and not using berendsens method...
   if( !( ( TemperatureReScaleTime < CONST_EPSILON ) && !( ConstantTemperature & 2 ) && ( Gamma >= CONST_EPSILON ) && ( ConstantTemperature & 16 ) ) )
   { 
    // calculate the temperature scaling factor...
    if( TheTemperature > 0.0 ) 
     ScaleFactor = sqrt( TargetTemperature / TheTemperature );
    else
     ScaleFactor = 1.0;
    
    // if needed and possible use Berendsen method 
    if( TemperatureReScaleTime >= CONST_EPSILON )         
     ScaleFactor = sqrt( 1.0 + MaxTimeStep * ( ScaleFactor * ScaleFactor - 1.0 ) / TemperatureReScaleTime );

    // re-initialize some data...
    TheAccumulatedTemperature -= TheTemperature;
    TheTotalKineticEnergy = 0.0;
    
    // now do the rescaling 
    for( int n = 0; n < Size; ++n )
     if( ( TheParticle = (*(TheState))[ n ] ) != NULL )
      if( TheParticle -> Mass > 0.0 )
      {
       TheParticle -> Momentum *= ScaleFactor;
       TheTotalKineticEnergy += 0.5 * Dot( TheParticle -> Momentum, TheParticle -> Momentum ) / TheParticle -> Mass;
      }

    // finaly re-calculate the temperature... 
    TheTemperature = ( BoltzmannConstant > 0.0 ? 0.66666666666666 * TheTotalKineticEnergy / ( theNrOfParticlesIncluded * BoltzmannConstant ) : 0.0 );
    TheAccumulatedTemperature += TheTemperature;
   }


  }
  
  return( 0 );
 }
 
    
 return( 1 );
}

/* ----------------------------------------------------------------------- */
/* Setup default values of the input data record                           */

void SetupGlobalInput( InputFileDataType& TheInputData )
{
 TheInputData.Particles = NULL;
 TheInputData.Interactions = NULL;
 TheInputData.Colors = NULL;
 TheInputData.Constraints = NULL;
 
 TheInputData.Time = 0.0;
 TheInputData.TimeStep = MAX_TIME_DELTA;
 TheInputData.EndTime = 0.0;
 TheInputData.AccuracyEnergy = ENERGY_CONSERVATION_CHECK;
 
 TheInputData.MaxStepSizeConformation = MAX_CONFORMATION_STEP_SIZE;
 TheInputData.AccuracyConformation = CONFORMATION_DEFAULT_ACCURACY;
 TheInputData.MaxStepsConformation = MAX_CONFORMATION_NSTEPS;
 
 TheInputData.DefaultParticle.Id = TheInputData.Periodicity = TheInputData.DoDynamics = TheInputData.DoConformation = 0;
 
 TheInputData.TheBox.X = TheInputData.TheBox.Y = TheInputData.TheBox.Z = 1.0;

 TheInputData.Scale.X = TheInputData.Scale.Y = TheInputData.Scale.Z = 1.0;
 TheInputData.Translation.X = TheInputData.Translation.Y = TheInputData.Translation.Z = 0.0;
 TheInputData.RotationAxis.X = TheInputData.RotationAxis.Y = TheInputData.RotationAxis.Z = 0.0;
 TheInputData.RotationAngle = 0.0;

 TheInputData.Periodicity = 0;
 TheInputData.TheBoxPointer = NULL;
 TheInputData.TargetTemperature = 0.0;
 TheInputData.GammaFactor = 0.0;
 TheInputData.ConstantTemperature = 0;
 TheInputData.BoltzmannConstant = 0.0;
 TheInputData.TemperatureReScaleTime = 0.0;
 TheInputData.RandomMomentaFactor = 0.0;
 TheInputData.CenterOfMassRemovalFactor = 0.0;
 
 TheInputData.NrOfAutoSnapshotsToTake = 0;
 TheInputData.AutoSnapShotTimeStep = MAX_TIME_DELTA;
 TheInputData.NextAutoSnapShotTime = MAX_TIME_DELTA;
}                            

/* ----------------------------------------------------------------------- */
/* Called by main program top set up the Dynamics object */

MyLangevinDynamicsPointer SetupDynamics( InputFileDataType& TheInputData ) 
{
 MyLangevinDynamicsPointer Dynamics;
 ParticlePointer         TheParticle;
 Vector                  CenterOfMassVelocity;
 double                  TotalMass, KineticEnergy, ScaleFactor;
 int                     NrOfParticles;

 if( ( fabs( TheInputData.Time ) < CONST_EPSILON ) && ( TheInputData.ConstantTemperature ) )
 {
  // init the velocities of the particles according to target temperature
  
  if( TheInputData.Particles != NULL )
  {
   KineticEnergy = TotalMass = CenterOfMassVelocity.X = CenterOfMassVelocity.Y = CenterOfMassVelocity.Z = 0.0;
   
   // first randomize all the momenta vectors
   for( int n = NrOfParticles = 0; n < TheInputData.Particles -> size(); ++n )
    if( ( TheParticle = (*(TheInputData.Particles))[ n ] ) != NULL )
    {
     TheParticle -> Momentum = TheParticle -> Mass * ARandomNormalizedVector();
     TotalMass += TheParticle -> Mass;     
     CenterOfMassVelocity += TheParticle -> Momentum;
     ++NrOfParticles;
    } 

    // calculate total center of mass velocity
    if( TotalMass > 0.0 )
     CenterOfMassVelocity /= TotalMass;
    else
     CenterOfMassVelocity *= 0.0;
    
    // remove total center of mass velocity from all particles to get total momentum 
    // to be zero and then calculate the current total kinetic energy of the system
    for( int n = 0; n < TheInputData.Particles -> size(); ++n )
     if( ( TheParticle = (*(TheInputData.Particles))[ n ] ) != NULL )
      if( TheParticle -> Mass > 0.0 )
      {
       TheParticle -> Momentum -= CenterOfMassVelocity * TheParticle -> Mass; 
       KineticEnergy += 0.5 * Dot( TheParticle -> Momentum, TheParticle -> Momentum ) / TheParticle -> Mass;
      } 

    // now determine scaling factor for velocities
    if( KineticEnergy > 0.0 ) 
     ScaleFactor = sqrt( ( 1.5 * NrOfParticles * TheInputData.TargetTemperature * TheInputData.BoltzmannConstant ) / KineticEnergy );
    else
     ScaleFactor = 1.0; 
      
    // now rescale the whole thing
    for( int n = 0; n < TheInputData.Particles -> size(); ++n )
     if( ( TheParticle = (*(TheInputData.Particles))[ n ] ) != NULL )
      TheParticle -> Momentum *= ScaleFactor;
           
  }
  
 }		
 
 Dynamics = new MyLangevinDynamics( TheInputData.Time, TheInputData.TimeStep, TheInputData.EndTime,
                                    TheInputData.TheBox, TheInputData.AccuracyEnergy, TheInputData.Periodicity,
		                    TheInputData.ConstantTemperature, TheInputData.TargetTemperature,
		                    TheInputData.BoltzmannConstant, TheInputData.TemperatureReScaleTime,
	    			    TheInputData.RandomMomentaFactor, TheInputData.GammaFactor, TheInputData.CenterOfMassRemovalFactor );
		  
				   
 if( Dynamics != NULL )
 {

  if( TheInputData.Particles != NULL )
   for( int n = 0; n < TheInputData.Particles -> size(); ++n )
    if( (*(TheInputData.Particles))[ n ] != NULL )
     Dynamics -> IncludeParticle( (*((*(TheInputData.Particles))[ n ])) );
     
 }
 
 return( Dynamics );
}

/* ----------------------------------------------------------------------- */
/* Called by main program top set up the steepestdescent conformation searcher object */

SteepestDescentPointer SetupConformation( InputFileDataType& TheInputData ) 
{
 SteepestDescentPointer Conformation;

 Conformation = new SteepestDescent();
 
 if( Conformation != NULL )
 {
 
  if( TheInputData.Particles != NULL )
   for( int n = 0; n < TheInputData.Particles -> size(); ++n )
    if( (*(TheInputData.Particles))[ n ] != NULL )
     Conformation -> IncludeParticle( (*((*(TheInputData.Particles))[ n ])) );
     
 }
 
 return( Conformation );
}

/* ----------------------------------------------------------------------- */
/* Stores the given state into the input data record                       */

void StoreFinalState( InputFileDataType &TheInputData, ParticleListPointer FinalState, int StoreMomentaToo = 0 )
{
 int n, Index;
    
 if( ( FinalState != NULL ) && ( TheInputData.Particles != NULL ) )
  for( n = 0; n < TheInputData.Particles -> size(); ++n )
   if( (*(TheInputData.Particles))[ n ] != NULL )
    for( Index = 0; Index < FinalState -> size(); ++Index )
     if( (*(FinalState))[ Index ] != NULL )
      if( (*(FinalState))[ Index ] -> Id == (*(TheInputData.Particles))[ n ] -> Id ) 
      {
       (*(TheInputData.Particles))[ n ] -> Position = (*(FinalState))[ Index ] -> Position;
       if( StoreMomentaToo ) (*(TheInputData.Particles))[ n ] -> Momentum = (*(FinalState))[ Index ] -> Momentum;
      } 

}

/* ----------------------------------------------------------------------- */
/* Gives the output to StdOutFilePointer and releases memory when run is   */
/* finished                                                                */

void DumpFinalStateAndCleanup( MyLangevinDynamicsPointer& Dynamics, InputFileDataType &TheInputData, int DoCleanUp = 1 )
{
 ParticleListPointer FinalState;
  
 // get the final state
 FinalState = Dynamics -> GetFinalStateParticlesPointer();
 
 if( ( FinalState != NULL ) && ( TheInputData.Particles != NULL ) )
 {
  if( DoCleanUp )
  { 
   fprintf( StdOutFilePointer, "\nDynamical simulation has finished:                     \n\n" );
   fprintf( StdOutFilePointer, "End time:                                              %.20le\n", Dynamics -> Time );
  }
  else
   fprintf( StdOutFilePointer, "\nTaking a snapshot:\n\nCurrent time:                                          %.20le\n", Dynamics -> Time ); 
   
  fprintf( StdOutFilePointer, "The current potential energy of system:                %.20le\n", Dynamics -> TheTotalInteractionPotential );
  fprintf( StdOutFilePointer, "The current kinetic energy of system:                  %.20le\n", Dynamics -> TheTotalKineticEnergy );
  fprintf( StdOutFilePointer, "The current temperature of the system:                 %.20le\n", Dynamics -> TheTemperature );
  fprintf( StdOutFilePointer, "The current pressure of the system:                    %.20le\n", Dynamics -> ThePressure );
  
  if( Dynamics -> TotalNrOfStepsDone >= 1 )
  {
   fprintf( StdOutFilePointer, "The average potential energy of system:                %.20le\n", ( Dynamics -> TheAccumulatedTotalInteractionPotential / Dynamics -> TotalNrOfStepsDone ) );
   fprintf( StdOutFilePointer, "The average kinetic energy of system:                  %.20le\n", ( Dynamics -> TheAccumulatedTotalKineticEnergy / Dynamics -> TotalNrOfStepsDone ) );  
   fprintf( StdOutFilePointer, "The average temperature of the system:                 %.20le\n", ( Dynamics -> TheAccumulatedTemperature / Dynamics -> TotalNrOfStepsDone ) );
   fprintf( StdOutFilePointer, "The average pressure of the system:                    %.20le\n", ( Dynamics -> TheAccumulatedPressure / Dynamics -> TotalNrOfStepsDone ) );
  }
  
  fprintf( StdOutFilePointer, "The density of system:                                 %.20le\n", ( TotalMass( (*FinalState) ) / ( Dynamics -> TheBox.X * Dynamics -> TheBox.Y * Dynamics -> TheBox.Z ) ) );
  fprintf( StdOutFilePointer, "Total mass of system:                                  %.20le\n", TotalMass( (*FinalState) ) );
  fprintf( StdOutFilePointer, "Total charge of system:                                %.20le\n", TotalCharge( (*FinalState) ) );
  fprintf( StdOutFilePointer, "Center of mass of system:                              [%.20le,%.20le,%.20le]\n", Dynamics -> TheCentreOfMass.X, Dynamics -> TheCentreOfMass.Y, Dynamics -> TheCentreOfMass.Z );
  fprintf( StdOutFilePointer, "Dipole moment of system:                               [%.20le,%.20le,%.20le]\n", Dynamics -> TheDipoleMoment.X, Dynamics -> TheDipoleMoment.Y, Dynamics -> TheDipoleMoment.Z );
  fprintf( StdOutFilePointer, "Mean of distance squared of system:                    %.20le\n\n", MeanDistanceSquared( (*FinalState) ) );
  
  // now copy the final state into the InputDataType for a possible restart dynamics run
  StoreFinalState( TheInputData, FinalState, 1 );
  TheInputData.TimeStep = Dynamics -> MaxTimeStep;
  TheInputData.Time = Dynamics -> Time;
  TheInputData.EndTime = Dynamics -> EndTime;
  TheInputData.TheBox = Dynamics -> TheBox;
  TheInputData.AccuracyEnergy = Dynamics -> Accuracy;
  TheInputData.Periodicity = Dynamics -> Periodicity;
  TheInputData.ConstantTemperature = Dynamics -> ConstantTemperature;
  TheInputData.BoltzmannConstant = Dynamics -> BoltzmannConstant;
  TheInputData.TargetTemperature = Dynamics -> TargetTemperature;
  TheInputData.TemperatureReScaleTime = Dynamics -> TemperatureReScaleTime;
  TheInputData.GammaFactor = Dynamics -> Gamma;
  TheInputData.RandomMomentaFactor = Dynamics -> RandomMomentaFactor;
  TheInputData.CenterOfMassRemovalFactor = Dynamics -> CenterOfMassRemovalFactor;

  // and dump the stuff into a file again
  if( !WriteDataToFile( OutputFileName, TheInputData, OutputFilePointer ) ) 
   fprintf( StdErrFilePointer, "Error writing to file '%s' !!\n\n", OutputFileName );
   
#if defined(USE_BOINC_API)
  CreateCheckPointFile();
#endif

  // remove the dynamics object
  if( DoCleanUp )
  {
#if !defined(USE_BOINC_GRAPH_API)
   delete Dynamics;     
   Dynamics = NULL;
  
   // and delete the data we wrote to the file
   DeleteTheData( TheInputData );
#endif
  } 
 }
  
}  

/* ----------------------------------------------------------------------- */
/* Saves the found constraint distances after minimizing to actual values  */

void SaveConstraintDistances( ConstraintListPointer Constraints, ParticleListPointer FinalState )
{
 DistanceConstraintPointer C;
 ParticlePointer P1, P2;
 int Index, P1ID, P2ID;
 
 if( ( Constraints == NULL ) || ( FinalState == NULL ) ) return;
 
 for( int n = 0; n < Constraints -> size(); ++n ) 
  if( ( C = (DistanceConstraintPointer)((*Constraints)[ n ]) ) != NULL )
   if( ( C -> Type == AGeneralDistanceConstraint ) && ( C -> Distance <= 0.0 ) )
   {
   
    // if it is a distance constraint we need to set, try to find the first two particles
    // in the constraint lista nd use those particles to determine the distance we need to set....
   
    for( Index = 0; ( Index < C -> GetParticleListSize() ) && ( !C -> GetParticleId( Index ) ); ++Index );
    P1ID = C -> GetParticleId( Index );
    
    for( Index++; ( Index < C -> GetParticleListSize() ) && ( !C -> GetParticleId( Index ) ); ++Index );
    P2ID = C -> GetParticleId( Index );
    
    if( P1ID & P2ID )
    {
     for( Index = 0; Index < FinalState -> size(); ++Index )
      if( ( P1 = (ParticlePointer)((*FinalState)[ Index ]) ) != NULL )
       if( P1 -> Id == P1ID ) break;
     
     for( Index = 0; Index < FinalState -> size(); ++Index )
      if( ( P2 = (ParticlePointer)((*FinalState)[ Index ]) ) != NULL )
       if( P2 -> Id == P2ID ) break;   
       
     C -> Distance = Abs( P1 -> Position - P2 -> Position );
    }
    else 
     C -> Distance = -1.0;
   }
  
}

/* ----------------------------------------------------------------------- */
/* Does the input reading and setting up of the dynamical run              */

int ReadInputDoConformationSearchAndSetupDynamics( void )
{
 ParticleListPointer FinalState;
 int Line = 0, N;
 
 // try and read input 
 if( !ReadDataFromFile( InputFileName, TheInputData, Line, InputFilePointer ) )
 {
  fprintf( StdErrFilePointer, "Error reading inputfile '%s' in line %d !\n", InputFileName, Line );
  DeleteTheData( TheInputData );
  return( 1 );
 }
 
 if( ( TheInputData.DoDynamics <= 0 && TheInputData.DoConformation <= 0 ) || TheInputData.Particles == NULL )
 {
  fprintf( StdErrFilePointer, "There is nothing to do for the program !\n\n" );
  DeleteTheData( TheInputData );
  return( 1 );
 }

 // we should do steepest descent before the possible dynamics run 
 if( TheInputData.DoConformation > 0 )
 {
 
  // setup conformation searcher class
  Conformation = SetupConformation( TheInputData );
  if( Conformation == NULL )
  {
   DeleteTheData( TheInputData );
   return( 1 );
  }

  fprintf( StdOutFilePointer, "Starting the conformational search:                    \n\n" );
  fprintf( StdOutFilePointer, "Max number of steps:                                   %d\n", TheInputData.MaxStepsConformation );
  fprintf( StdOutFilePointer, "Max step size:                                         %.20le\n", TheInputData.MaxStepSizeConformation );
  fprintf( StdOutFilePointer, "Max error for convergence:                             %.20le\n\n", TheInputData.AccuracyConformation );
   
  fprintf( StdOutFilePointer, "Total potential energy of system at start of search:   %.20le\n", TotalPotentialEnergy( (*(Conformation -> GetFinalStateInteractionsPointer())) ) );
  fprintf( StdOutFilePointer, "Mean of distance squared of system at start of search: %.20le\n\n", MeanDistanceSquared( (*(Conformation -> GetFinalStateParticlesPointer())) ) );
  
  N = Conformation -> FindMinima( TheInputData.MaxStepsConformation, TheInputData.MaxStepSizeConformation, TheInputData.AccuracyConformation, TheInputData.Constraints );

  FinalState = Conformation -> GetFinalStateParticlesPointer();
    
  fprintf( StdOutFilePointer, "Conformation searcher stopped at itteration:           %d\n\n", N );
  fprintf( StdOutFilePointer, "Total potential energy of system:                      %.20le\n\n",TotalPotentialEnergy( (*(Conformation -> GetFinalStateInteractionsPointer())) ) );
  fprintf( StdOutFilePointer, "Center of mass of system:                              [%.20le,%.20le,%.20le]\n",   CentreOfMass( (*FinalState) ).X, CentreOfMass( (*FinalState) ).Y, CentreOfMass( (*FinalState) ).Z );
  fprintf( StdOutFilePointer, "Dipole moment of system:                               [%.20le,%.20le,%.20le]\n", DipoleMoment( (*FinalState) ).X, DipoleMoment( (*FinalState) ).Y, DipoleMoment( (*FinalState) ).Z );
  fprintf( StdOutFilePointer, "Mean of distance squared of system:                    %.20le\n\n", MeanDistanceSquared( (*FinalState) ) );
  
  // now copy the found conformation into the InputDataType for a possible dynamics run
  FinalState = Conformation -> GetFinalStateParticlesPointer();
  StoreFinalState( TheInputData, FinalState, 0 );
   
  // Also set all distance constraints to their actual values we found by minimizing...
  SaveConstraintDistances( TheInputData.Constraints, FinalState );
  
  TheInputData.DoConformation = -1;
       
  // and dump the stuff into a file again
  if( !WriteDataToFile( OutputFileName, TheInputData, OutputFilePointer ) ) 
   fprintf( StdErrFilePointer, "Error writing to file '%s'!!\n\n", OutputFileName );

  // and finally remove conformation searcher object 
  delete Conformation;
  Conformation = NULL;
 }
 
 
 // we should do a dynamics run, so set it up
 if( TheInputData.DoDynamics > 0 )
 {
 
  // setup the dynamics class  
  Dynamics = SetupDynamics( TheInputData );
  if( Dynamics == NULL )
  {
   DeleteTheData( TheInputData );
   return( 1 );
  }

  fprintf( StdOutFilePointer, "Starting the dynamical simulation:                     \n\n" );
  if( Dynamics -> Periodicity )
   fprintf( StdOutFilePointer, "Periodic box:                                          [%.20le,%.20le,%.20le]\n", Dynamics -> TheBox.X, Dynamics -> TheBox.Y, Dynamics -> TheBox.Z );
  else
   fprintf( StdOutFilePointer, "Single cell:                                           [%.20le,%.20le,%.20le]\n", Dynamics -> TheBox.X, Dynamics -> TheBox.Y, Dynamics -> TheBox.Z );

  fprintf( StdOutFilePointer, "Start time:                                            %.20le\n", Dynamics -> StartTime );
  fprintf( StdOutFilePointer, "End time:                                              %.20le\n", Dynamics -> EndTime );
  fprintf( StdOutFilePointer, "Max time step:                                         %.20le\n", Dynamics -> MaxTimeStep );
  fprintf( StdOutFilePointer, "Number of automatic snapshots to take in this run:     %u\n", TheInputData.NrOfAutoSnapshotsToTake );
  
  if( Dynamics -> ConstantTemperature <= 0 )
   fprintf( StdOutFilePointer, "Max error in conservation of total energy:             %.20le\n", Dynamics -> Accuracy );
   
  if( Dynamics -> ConstantTemperature )
  {
   fprintf( StdOutFilePointer, "The boltzmann constant:                                %.20le\n", Dynamics -> BoltzmannConstant );
   if( Dynamics -> ConstantTemperature < 0 && ( fabs( Dynamics -> Time ) < CONST_EPSILON ) )
    fprintf( StdOutFilePointer, "The temperature will be initialized to:                %.20le\n", Dynamics -> TargetTemperature );
   if( Dynamics -> ConstantTemperature > 0 )
   {
    fprintf( StdOutFilePointer, "The temperature will be kept constant at:              %.20le\n", Dynamics -> TargetTemperature );
    fprintf( StdOutFilePointer, "The temperature rescaling time:                        %.20le\n", Dynamics -> TemperatureReScaleTime );
    fprintf( StdOutFilePointer, "The temperature random momenta factor:                 %.20le\n", Dynamics -> RandomMomentaFactor );     
    fprintf( StdOutFilePointer, "The temperature center of mass removal factor:         %.20le\n", Dynamics -> CenterOfMassRemovalFactor );     
    fprintf( StdOutFilePointer, "The temperature langevin gamma factor:                 %.20le\n", Dynamics -> Gamma );     
   } 
  } 
 }

 return( 0 );
}

/* ----------------------------------------------------------------------- */
/* Redirects StdErrFilePointer and StdOutFilePointer file pointers if needed ...                 */

void RedirectStdOutAndStdErr( void )
{
 if( StdOutFileName[ 0 ]  )
 {
  StdOutFilePointer = fopen( StdOutFileName, "wt" );
  if( StdOutFilePointer == NULL )
  {
   fprintf( stderr, "Error opening stdout filename !\n" );
   exit( 1 );
  }
 } 
 else
  StdOutFilePointer = stdout;


 if( StdErrFileName[ 0 ]  )
 {
  StdErrFilePointer = fopen( StdErrFileName, "wt" );
  if( StdErrFilePointer == NULL )
  {
   fprintf( stderr, "Error opening stderr filename !\n" );
   exit( 1 );
  }
 }
 else
  StdErrFilePointer = stderr;

}

/* ----------------------------------------------------------------------- */
/* Does a a few timesteps of work                                          */

void DoSomeOfTheWork( MyLangevinDynamicsPointer Dynamics )
{

  if( ( Dynamics -> EndTime - Dynamics -> Time ) < CONST_EPSILON * Dynamics -> EndTime )
  {
   Dynamics -> StopTrajectory = 1;
   return;
  } 
    
  if( ( Dynamics -> MaxTimeStep - ( Dynamics -> Time - Dynamics -> PreviousTime ) ) < CONST_EPSILON * Dynamics -> MaxTimeStep )
  {
   Dynamics -> NrOfTimeStepsDone = 0;
   Dynamics -> PreviousTime = Dynamics -> Time;
   Dynamics -> TimeStep = Dynamics -> MaxTimeStep;
  }

  // only find appropriate integration timesteps for conservation of energy 
  // when we are not rescaling velocities for constant T dynamics ...
  if( Dynamics -> ConstantTemperature <= 0 )
  {
   Dynamics -> TimeStep = Dynamics -> NewTimeStep;
   Dynamics -> NrOfStepsToDo = ( int )( floor( 0.25 * ( MIN( Dynamics -> MaxTimeStep, Dynamics -> EndTime - Dynamics -> Time ) - ( Dynamics -> Time - Dynamics -> PreviousTime ) ) / Dynamics -> TimeStep ) );
   Dynamics -> NrOfStepsToDo = ( Dynamics -> NrOfStepsToDo >= 1 ? Dynamics -> NrOfStepsToDo : 1 );
   Dynamics -> FindAppropriateTimeStep = 1;  
  }
  else
   Dynamics -> NrOfStepsToDo = 1;
   
  // if we do take automatic snapshots... 
  if( TheInputData.NrOfAutoSnapshotsToTake > 0 ) 
   if( ( TheInputData.NextAutoSnapShotTime - Dynamics -> Time ) < CONST_EPSILON * TheInputData.NextAutoSnapShotTime )
   {
    TheInputData.NextAutoSnapShotTime = TheInputData.NextAutoSnapShotTime + TheInputData.AutoSnapShotTimeStep;
    --TheInputData.NrOfAutoSnapshotsToTake;

    // take snapshot now ...
    DumpFinalStateAndCleanup( Dynamics, TheInputData, 0 );
   }
  
  
  Dynamics -> NrOfTimeStepsDone += Dynamics -> RunTrajectory( Dynamics -> Time, Dynamics -> TimeStep, Dynamics -> NrOfStepsToDo, TheInputData.Constraints );

}  

/* ----------------------------------------------------------------------- */
/* Stuff to get GLUT / BOINC working                                       */

#if defined(USE_GLUT_GRAPH_API) || defined(USE_BOINC_GRAPH_API)

#   include "OpenGl.inc"

#   if defined(USE_GLUT_GRAPH_API)
#      include "Glut.inc"
#   elif defined(USE_BOINC_GRAPH_API)
#      include "Boinc.inc"
#    endif

#elif defined(USE_BOINC_API)

#   include "Boinc.inc"

#endif

/* ------------------------------------------------------- Main program -- */

int main(  int argc, char *argv[] ) 
{
 int Error;

#if defined(USE_CRLIBM)
 crlibm_init();              // init the crlibm library 
#endif
 
 // set random number seed to current time in seconds
 SRand( time( NULL ) );

// setup default values for the data
 SetupGlobalInput( TheInputData );

 if( argc >= 2 )
  strcpy( InputFileName, argv[ 1 ] );
 else
  strcpy( InputFileName, "classical.in" );
  
 if( argc >= 3 )
  strcpy( OutputFileName, argv[ 2 ] ); 
 else 
  strcpy( OutputFileName, "classical.out" ); 
  
 if( argc >= 4 )
  strcpy( StdOutFileName, argv[ 3 ] ); 
 else 
  StdOutFileName[ 0 ] = '\0';

 if( argc >= 5 )
  strcpy( StdErrFileName, argv[ 4 ] ); 
 else 
  StdErrFileName[ 0 ] = '\0';

#if defined(USE_GLUT_GRAPH_API)

  // setup StdOutFilePointer and StdErrFilePointer...
  RedirectStdOutAndStdErr();

  // the setup the glut context and start rendering the crap
  Error = ReadInputDoConformationSearchAndSetupDynamics();
  
  if( Error ) return( Error );
  
  if( TheInputData.DoDynamics )
   SetupGlutRendering( argc, argv );

  if( StdOutFileName[ 0 ] ) fclose( StdOutFilePointer );   // only close if they have been redirected ...
  if( StdErrFileName[ 0 ] ) fclose( StdErrFilePointer );
  
#elif defined(USE_BOINC_GRAPH_API)

  // run graphical boinc  
  SetupBoincRendering( argc, argv );
  
#elif defined(USE_BOINC_API)

  // run non-graphical boinc
  SetupBoinc( argc, argv );  

#elif defined(INPUT_FILE_CHECKER)

  int Line = 0;

  // setup StdOutFilePointer and StdErrFilePointer...
  RedirectStdOutAndStdErr();

  // try and read input
  Error = !ReadDataFromFile( InputFileName, TheInputData, Line, InputFilePointer );

  if( Error )
   fprintf( StdErrFilePointer, "Error reading inputfile '%s' in line %d !\n", InputFileName, Line );
  
  DeleteTheData( TheInputData );

  if( StdOutFileName[ 0 ] ) fclose( StdOutFilePointer );   // only close if they have been redirected ...
  if( StdErrFileName[ 0 ] ) fclose( StdErrFilePointer );

  if( Error )
   return( Line );

#else  
  
  // setup StdOutFilePointer and StdErrFilePointer...
  RedirectStdOutAndStdErr();

  // start running the trajectory, quietly 
  Error = ReadInputDoConformationSearchAndSetupDynamics();
  
  if( Error ) return( Error );
  
  if( TheInputData.DoDynamics )
  {
   while( !Dynamics -> StopTrajectory )
    DoSomeOfTheWork( Dynamics );
    
   DumpFinalStateAndCleanup( Dynamics, TheInputData );
  }

  if( StdOutFileName[ 0 ] ) fclose( StdOutFilePointer );   // only close if they have been redirected ...
  if( StdErrFileName[ 0 ] ) fclose( StdErrFilePointer );
  
#endif
 
 return( 0 );
}