/*   []----------------------------------------------[]
      |              ClassicalDynamics.h             |
     []----------------------------------------------[]
      |                                              |
      |   AUTHOR:    MFSomers 2005.                  |
      |   USE:       Defines a classical dynamics    |
      |              framework class to be used for  |
      |              different equations of motions  |
      |              integrators... It also defines  |
      |              a particle and an interaction   |
      |              class to be used in this        |
      |              framework...                    |
      |                                              |
     []----------------------------------------------[]
*/

// 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 Lesser General Public License as published by the Free Software Foundation.
//
// http://www.gnu.org/licenses/lgpl.txt

#if !defined(__CLASSICAL_DYNAMICS_INCLUDED__)

#define __CLASSICAL_DYNAMICS_INCLUDED__

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

#include <vector>

#define DEFAULT_NR_OF_STATES          2

#define DEFAULT_START_NR_INTERACTIONS 32
#define DEFAULT_START_NR_PARTICLES    128 

#define DEFAULT_DELTA_NR_INTERACTIONS 16
#define DEFAULT_DELTA_NR_PARTICLES    64

/* ----------------------------------------------------------------------- */
/* Define forwards the base classes defined in this header */

class Particle;
class Interaction;
class Constraint;
class ClassicalStates;
class ClassicalDynamics;

typedef class Particle    *ParticlePointer;
typedef class Interaction *InteractionPointer;
typedef class Constraint  *ConstraintPointer;
typedef class ClassicalStates *ClassicalStatesPointer;
typedef class ClassicalDynamics *ClassicalDynamicsPointer;

typedef vector<ParticlePointer> ParticleList;
typedef vector<InteractionPointer> InteractionList;
typedef vector<ConstraintPointer> ConstraintList;

typedef ParticleList *ParticleListPointer;
typedef InteractionList *InteractionListPointer;
typedef ConstraintList *ConstraintListPointer;

typedef vector<InteractionListPointer> ListOfInteractionLists;
typedef vector<ParticleListPointer> ListOfParticleLists;

/* ----------------------------------------------------------------------- */
/* Definition of a particle */

class Particle { friend class ClassicalStates;
                 public:

                 Particle( unsigned int TheID = GET_NEW_ID );
                 Particle( Vector ThePosition, Vector TheMomentum, double TheMass = 0.0, double TheCharge = 0.0,
                           double TheRadius = 0.0, int TheColorIndex = 0, unsigned int TheID = GET_NEW_ID );
		 Particle( Particle& P );
		 virtual ~Particle()
                 {
		  ExcludeAllInteractions();
                 }
		 
		 virtual Particle& operator=( Particle& P );
		 virtual ParticlePointer Clone()
		 {
		  return( new Particle( (*this) ) );
		 }
		 virtual ParticlePointer Create( unsigned int TheID = GET_NEW_ID )
		 {
		  return( new Particle( TheID ) );
		 }
		 
		 void IncludeInteraction( Interaction& I );
		 void ExcludeInteraction( unsigned int ID );
		 void ExcludeInteraction( Interaction& I );
		 void DeleteInteraction( unsigned int ID );
		 void DeleteInteraction( Interaction& I );		 
		 int  IsInteractionIncluded( unsigned int ID );
		 int  IsInteractionIncluded( Interaction& I );
		 InteractionPointer GetInteractionPointer( int Index );
		 void ExcludeAllInteractions( void );
		 void DeleteAllInteractions( void );
		 
		 double TheTotalInteractionPotential( void );
		 Vector TheTotalForce( void ); 
		 
		 int NrOfInteractionsIncluded( void );
		 int GetInteractionListSize( void );
		 
		 inline double TheKineticEnergy( void )
		 {
		  return( Mass != 0.0 ? 0.5 * Dot( Momentum, Momentum ) / Mass : 0.0 ); 
		 } 	
		 
		 inline double TheTotalEnergy( void )
		 {
		  return( TheKineticEnergy() + TheTotalInteractionPotential() );
		 } 
		 
                 // now the data members...
		 
                 Vector                 Position;
                 Vector                 Momentum;
                 double                 Mass;
                 double                 Charge;
                 double                 Radius;
		 unsigned int           Id;
		 int                    TheInternalFlag;
		 int                    ColorIndex;
		 ParticleType           Type;
		 		 
		 protected:
		 
		 InteractionList        Interactions;
		 ClassicalStatesPointer DynamicsContext;
               };
	       
/* ----------------------------------------------------------------------- */   
/* Definition of a general interaction */
       
class Interaction { friend class ClassicalStates;
                    public:
		    
		    Interaction( unsigned int TheId = GET_NEW_ID, Vector *UsePeriocBox = NULL );
	            Interaction( Interaction& I );
		    virtual ~Interaction()
		    {
		     ExcludeAllParticles();
		    }
		    
		    virtual Interaction& operator=( Interaction& I );
		    virtual InteractionPointer Clone()
		    {
		     return( new Interaction( (*this) ) );
		    }
		    virtual InteractionPointer Create( unsigned int TheId = GET_NEW_ID, Vector *UsePeriocBox = NULL )
		    {
		     return( new Interaction( TheId, UsePeriocBox ) );
		    }
		    
		    void IncludeParticle( Particle& P );
	            void ExcludeParticle( unsigned int ID );
	            void ExcludeParticle( Particle& P );
		    void DeleteParticle( unsigned int ID );
		    void DeleteParticle( Particle& P );		    
		    int  IsParticleIncluded( unsigned int ID );
		    int  IsParticleIncluded( Particle& P );
		    ParticlePointer GetParticlePointer( int Index );
		    void ExcludeAllParticles( void );
		    void DeleteAllParticles( void );
                 
		    int NrOfParticlesIncluded( void );
		    int GetParticleListSize( void );
		    
		    virtual double InteractionPotential( void ); 
		    virtual Vector InteractionForce( unsigned int ID );    
		    Vector InteractionForce( Particle &P );
	    
                    unsigned int           Id;
		    int                    TheInternalFlag;
        	    InteractionType        Type;
		    Vector                 *ThePeriodicBoxPointer;
		    		   		    		    
		    protected:
		    
		    ParticleList           Particles;
		    ClassicalStatesPointer DynamicsContext;
                  };

/* ----------------------------------------------------------------------- */
/* Definition of a general constraint */
		  		  
class Constraint { public:

                   Constraint( unsigned int TheID = GET_NEW_ID );
		   Constraint( Constraint& C );
		   
		   virtual Constraint& operator=( Constraint& C );
		   
		   void IncludeParticle( Particle& P );
	           void ExcludeParticle( unsigned int ID );
	           void ExcludeParticle( Particle& P );
	           int  IsParticleIncluded( unsigned int ID );
		   int  IsParticleIncluded( Particle& P );
		   unsigned int GetParticleId( int Index );
 	           void ExcludeAllParticles( void );
             
	           int NrOfParticlesIncluded( void );
		   int GetParticleListSize( void );
		   
		   unsigned int         Id;
		   ConstraintType       Type;
                   int                  TheInternalFlag;		   

		   vector<unsigned int> ParticleIDs;
		   vector<int>          ParticleIndices;
		 };

/* ----------------------------------------------------------------------- */
/* Definition of a general collection of classical states */

class ClassicalStates { public:
			  
			ClassicalStates( int NrOfStates = DEFAULT_NR_OF_STATES, unsigned int TheId = GET_NEW_ID );
			ClassicalStates( ClassicalStates& D );
			virtual ~ClassicalStates()
			{
			 DeleteAllTheParticleLists();
			 DeleteAllTheInteractionLists();
			}
			  
			virtual ClassicalStates& operator=( ClassicalStates& D );
			virtual ClassicalStatesPointer Clone()
			{
			 return( new ClassicalStates( (*this) ) );
			}
			virtual ClassicalStatesPointer Create( int NrOfStates = DEFAULT_NR_OF_STATES, unsigned int TheId = GET_NEW_ID  )
			{
			 return( new ClassicalStates( NrOfStates, TheId ) );
			}			  
			  
			void IncludeParticle( Particle& P );
			void IncludeInteraction( Interaction& I );

			int  NrOfInteractionsIncluded( void );
			int  NrOfParticlesIncluded( void );
			int  NrOfStatesIncluded( void );
			 
		        int  IsInteractionIncluded( unsigned int ID );
		        int  IsInteractionIncluded( Interaction& I );			  
		        int  IsParticleIncluded( unsigned int ID );
		        int  IsParticleIncluded( Particle& P );			  
			  			  
			void DeleteParticle( Particle& P );
			void DeleteInteraction( Interaction& I );
			void DeleteParticle( unsigned int ID );
			void DeleteInteraction( unsigned int ID );
			  
			void DeleteAllParticles( void );
			void DeleteAllInteractions( void );
			  
			void IncludeParticleInInteraction( Particle& P, Interaction& I );
			void IncludeInteractionInParticle( Interaction& I, Particle& P );			  
			  
			void ExcludeInteractionFromParticle( unsigned int IID, unsigned int PID );
			void ExcludeParticleFromInteraction( unsigned int PID, unsigned int IID );
			void ExcludeInteractionFromParticle( Interaction& I, Particle& P );
			void ExcludeParticleFromInteraction( Particle& P, Interaction& I );			  

			void ExcludeAllInteractionsFromParticle( unsigned int PID );
			void ExcludeAllParticlesFromInteraction( unsigned int IID );
			void ExcludeAllInteractionsFromParticle( Particle& P );
			void ExcludeAllParticlesFromInteraction( Interaction& I );			  
			  
			void DeleteAllInteractionsOfParticle( unsigned int PID );
			void DeleteAllParticlesOfInteraction( unsigned int IID );
			void DeleteAllInteractionsOfParticle( Particle& P );
			void DeleteAllParticlesOfInteraction( Interaction& I );			  		  
			  
			ParticleListPointer    GetStateParticlesPointer( int s );
			InteractionListPointer GetStateInteractionsPointer( int s );
			ParticlePointer     GetParticlePointer( int Index, int s );
			InteractionPointer  GetInteractionPointer( int Index, int s );
			int                 GetParticleListSize( void );
			int                 GetInteractionListSize( void );
			int                 GetStateListSize( void );
			  
			void CycleTheStates( int Steps = -1 );
			  
			// the data members 
			 
			unsigned int           Id;
			int                    TheInternalFlag;
			ClassicalDynamicsType  Type;
			  			  
			protected:

			void DeleteAllTheParticleLists( void );
			void DeleteAllTheInteractionLists( void );				  
			  			  
			ListOfParticleLists    TheParticleListList;
			ListOfInteractionLists TheInteractionListList;
                      };
			   
/* ----------------------------------------------------------------------- */
/* A general ClassicalDynamics class definition */

class ClassicalDynamics : public ClassicalStates
                        { public:
			  
			  ClassicalDynamics( int NrOfStates = DEFAULT_NR_OF_STATES, unsigned int TheId = GET_NEW_ID ) : 
			                     ClassicalStates( NrOfStates, TheId )
			  {
			   Type = GeneralClassicalDynamics;
			  };
			  ClassicalDynamics( ClassicalDynamics& D ) : ClassicalStates( D )
			  {
			  };
			  
			  virtual ClassicalDynamics& operator=( ClassicalDynamics& D )
			  {
			   if( &D == &(*this) ) return( (*this) );
			   ClassicalStates::operator=( D );
			   return( (*this) );
			  };
			  
			  virtual ClassicalDynamicsPointer Clone()
			  {
			   return( new ClassicalDynamics( (*this) ) );
			  }
			  virtual ClassicalDynamicsPointer Create( int NrOfStates = DEFAULT_NR_OF_STATES, unsigned int TheId = GET_NEW_ID  )
			  {
			   return( new ClassicalDynamics( NrOfStates, TheId ) );
			  }			  
			  
			  virtual ParticleListPointer GetFinalStateParticlesPointer( void )
			  {
			   return( GetStateParticlesPointer( 1 ) );
			  }
			  
			  virtual InteractionListPointer GetFinalStateInteractionsPointer( void )
			  {
			   return( GetStateInteractionsPointer( 1 ) );
			  }			  
			  
			  virtual int Integrator( double Step, double Start = 0.0, ConstraintListPointer Constraints = NULL )
			  {
			   return( 0 );
			  } 
			  
			  virtual int IntegratorConservationCheck( double Step )
			  {
			   return( 0 );
			  } 
			  
			  virtual int TrajectoryContinuationCheck( double At, double Step, int nStep, int nSteps )
			  {
			   return( 0 );
			  } 
			  
			  virtual int RunTrajectory( double Start, double Step, int nSteps, ConstraintListPointer Constraints = NULL );
	                };
			  
/* ----------------------------------------------------------------------- */

Vector CentreOfMass( const ParticleList& TheParticles );
Vector DipoleMoment( const ParticleList& TheParticles );
double TotalPotentialEnergy( const InteractionList& TheInteractions );
double TotalKineticEnergy( const ParticleList& TheParticles );
double TotalMass( const ParticleList& TheParticles );
double TotalCharge( const ParticleList& TheParticles );
double Virial( const ParticleList& TheParticles, const Vector& TheBox );
double MeanDistanceSquared( const ParticleList& TheParticles );

#endif