/*   []----------------------------------------[]
      |                FiniteDiff.h            |
     []----------------------------------------[]
      |                                        |
      |   AUTHOR:    MFSomers 2005.            |
      |   USE:       Finite differences on     |
      |               interactions...          |      
      |                                        |
     []----------------------------------------[]
*/

// 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(__FINITEDIFF_INCLUDED__)

#define __FINITEDIFF_INCLUDED__

#include "Global.h"
#include "Vector.h"
#include "Constants.h"
#include "ClassicalDynamics.h"
#include "PointPair.h"
#include "Bending.h"
#include "Torsional.h"

/* ----------------------------------------------------------------------- */
/* Definition of an interaction that calculates forces using finite differences */

class   FiniteDifferenceInteraction;
typedef FiniteDifferenceInteraction *FiniteDifferenceInteractionPointer;

class FiniteDifferenceInteraction : public Interaction
                                    {
				     public:
				    
				     inline FiniteDifferenceInteraction( double TheDelta = CONST_EPSILON, 
				              unsigned int TheId = GET_NEW_ID, Vector *UsePeriodicBox = NULL ): Interaction( TheId, UsePeriodicBox )
				     {
				      DifferenceDelta = ( TheDelta < CONST_EPSILON ? CONST_EPSILON : TheDelta );
				      Type = AGeneralFiniteDifferenceInteraction; 
				     }
				     inline FiniteDifferenceInteraction( FiniteDifferenceInteraction& I ) : Interaction( I )
				     {
				      DifferenceDelta = I.DifferenceDelta;
				     }

				     virtual FiniteDifferenceInteraction& operator=( FiniteDifferenceInteraction& I )
				     {
				      if( &I == &(*this) ) return( (*this) );
				      Interaction::operator=( I );
				      DifferenceDelta = I.DifferenceDelta;
				      return( (*this) );
				     }				     
		                     virtual FiniteDifferenceInteractionPointer Clone()
		                     {
		                      return( new FiniteDifferenceInteraction( (*this) ) );
		                     }
		                     virtual FiniteDifferenceInteractionPointer Create(double TheDelta = CONST_EPSILON, unsigned int TheId = GET_NEW_ID, Vector *UsePeriodicBox = NULL )
		                     {
		                      return( new FiniteDifferenceInteraction( TheDelta, TheId, UsePeriodicBox ) );
		                     }				     
				     
				     virtual Vector InteractionForce( unsigned int ID ); 
				     
				     protected:
				     
				     double DifferenceDelta;
                                    }; 

/* ----------------------------------------------------------------------- */    
/* A point pair interaction that calculates the force through finite differences */

class   FiniteDifferencePointPairInteraction;
typedef FiniteDifferencePointPairInteraction *FiniteDifferencePointPairInteractionPointer;

class FiniteDifferencePointPairInteraction : public PointPairInteraction
                                             {
					      public:
					      
				              inline FiniteDifferencePointPairInteraction( double TheDelta = CONST_EPSILON, 
				                     unsigned int TheId = GET_NEW_ID, Vector *UsePeriodicBox = NULL ): PointPairInteraction( TheId, UsePeriodicBox )
					      {
				               DifferenceDelta = ( TheDelta < CONST_EPSILON ? CONST_EPSILON : TheDelta );
					       Type = AGeneralFiniteDifferencePointPairInteraction;
					      }
					      inline FiniteDifferencePointPairInteraction( FiniteDifferencePointPairInteraction& I ) : PointPairInteraction( I )
					      {
					       DifferenceDelta = I.DifferenceDelta;
					      }

					      virtual FiniteDifferencePointPairInteraction& operator=( FiniteDifferencePointPairInteraction& I )
					      {
					       PointPairInteraction::operator=( I );
					       DifferenceDelta = I.DifferenceDelta;
					       return( (*this) );
					      } 
		                              virtual FiniteDifferencePointPairInteractionPointer Clone()
		                              {
		                               return( new FiniteDifferencePointPairInteraction( (*this) ) );
		                              }
		                              virtual FiniteDifferencePointPairInteractionPointer Create( double TheDelta = CONST_EPSILON, 
					                               unsigned int TheId = GET_NEW_ID, Vector *UsePeriodicBox = NULL )
		                              {
		                               return( new FiniteDifferencePointPairInteraction( TheDelta, TheId, UsePeriodicBox ) );
		                              }					      
					      
					      virtual Vector PointPairForce( Particle& P1, Particle& P2 );
					      
					      protected:
					      
					      double DifferenceDelta;
					     };

/* ----------------------------------------------------------------------- */    
/* A bond bending interaction that calculates the force through finite differences */

class   FiniteDifferenceBendingInteraction;
typedef FiniteDifferenceBendingInteraction *FiniteDifferenceBendingInteractionPointer;

class FiniteDifferenceBendingInteraction : public BendingInteraction
                                             {
					      public:
					      
				              inline FiniteDifferenceBendingInteraction( double TheDelta = CONST_EPSILON, 
				                     unsigned int TheId = GET_NEW_ID, Vector *UsePeriodicBox = NULL ): BendingInteraction( TheId, UsePeriodicBox )
					      {
				               DifferenceDelta = ( TheDelta < CONST_EPSILON ? CONST_EPSILON : TheDelta );
					       Type = AGeneralFiniteDifferenceBendingInteraction;
					      }
					      inline FiniteDifferenceBendingInteraction( FiniteDifferenceBendingInteraction& I ) : BendingInteraction( I )
					      {
					       DifferenceDelta = I.DifferenceDelta;
					      }

					      virtual FiniteDifferenceBendingInteraction& operator=( FiniteDifferenceBendingInteraction& I )
					      {
					       BendingInteraction::operator=( I );
					       DifferenceDelta = I.DifferenceDelta;
					       return( (*this) );
					      } 
		                              virtual FiniteDifferenceBendingInteractionPointer Clone()
		                              {
		                               return( new FiniteDifferenceBendingInteraction( (*this) ) );
		                              }
		                              virtual FiniteDifferenceBendingInteractionPointer Create( double TheDelta = CONST_EPSILON, 
				                     unsigned int TheId = GET_NEW_ID, Vector *UsePeriodicBox = NULL )
		                              {
		                               return( new FiniteDifferenceBendingInteraction( TheDelta, TheId, UsePeriodicBox ) );
		                              }					      
					      
					      virtual double DerivativeOfBendingPotential( double cosangle );
					      					      
					      protected:
					      
					      double DifferenceDelta;
					     };
					     
/* ----------------------------------------------------------------------- */    
/* A bond torsional interaction that calculates the force through finite differences */

class   FiniteDifferenceTorsionalInteraction;
typedef FiniteDifferenceTorsionalInteraction *FiniteDifferenceTorsionalInteractionPointer;

class FiniteDifferenceTorsionalInteraction : public TorsionalInteraction
                                             {
					      public:
					      
				              inline FiniteDifferenceTorsionalInteraction( double TheDelta = CONST_EPSILON, 
				                     unsigned int TheId = GET_NEW_ID, Vector *UsePeriodicBox = NULL ): TorsionalInteraction( TheId, UsePeriodicBox )
					      {
				               DifferenceDelta = ( TheDelta < CONST_EPSILON ? CONST_EPSILON : TheDelta );
					       Type = AGeneralFiniteDifferenceTorsionalInteraction;
					      }
					      inline FiniteDifferenceTorsionalInteraction( FiniteDifferenceTorsionalInteraction& I ) : TorsionalInteraction( I )
					      {
					       DifferenceDelta = I.DifferenceDelta;
					      }

					      virtual FiniteDifferenceTorsionalInteraction& operator=( FiniteDifferenceTorsionalInteraction& I )
					      {
					       TorsionalInteraction::operator=( I );
					       DifferenceDelta = I.DifferenceDelta;
					       return( (*this) );
					      } 
		                              virtual FiniteDifferenceTorsionalInteractionPointer Clone()
		                              {
		                               return( new FiniteDifferenceTorsionalInteraction( (*this) ) );
		                              }
		                              virtual FiniteDifferenceTorsionalInteractionPointer Create( double TheDelta = CONST_EPSILON, 
				                         unsigned int TheId = GET_NEW_ID, Vector *UsePeriodicBox = NULL )
		                              {
		                               return( new FiniteDifferenceTorsionalInteraction( TheDelta, TheId, UsePeriodicBox ) );
		                              }					      
					      
					      virtual double DerivativeOfTorsionalPotential( double cosangle );
					      					      
					      protected:
					      
					      double DifferenceDelta;
					     };
					     
/* ----------------------------------------------------------------------- */    

#endif