/*   []----------------------------------------[]
      |             FiniteDiff.cpp             |
     []----------------------------------------[]
      |                                        |
      |   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

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

#include "FiniteDiff.h"

/* ----------------------------------------------------------------------- */

Vector FiniteDifferenceInteraction::InteractionForce( unsigned int ID )
{
 ParticlePointer   TheParticle;
 Vector            TheForce, TheOldPosition;
 int               ParticleNumber;
 double            ScaledStep;
 
 TheForce.X = TheForce.Y = TheForce.Z = 0.0;
 
 // If particle is not feeling this interaction, return zero force...
 if( ( ParticleNumber = IsParticleIncluded( ID ) ) <= -1 )
  return( TheForce );
  
 TheParticle = Particles[ ParticleNumber ];
 TheOldPosition = TheParticle -> Position; 
 
 // calculate derivative along X through threepoint finite difference... 
 ScaledStep = fabs( TheOldPosition.X ) * DifferenceDelta;
 TheParticle -> Position.X += ScaledStep;
 TheForce.X = InteractionPotential();
 TheParticle -> Position.X -= ScaledStep;
 TheParticle -> Position.X -= ScaledStep;
 TheForce.X = 0.5 * ( InteractionPotential() - TheForce.X ) / ScaledStep; 
 TheParticle -> Position.X += ScaledStep;
 
 // calculate derivative along Y through threepoint finite difference... 
 ScaledStep = fabs( TheOldPosition.Y ) * DifferenceDelta;
 TheParticle -> Position.Y += ScaledStep;
 TheForce.Y = InteractionPotential();
 TheParticle -> Position.Y -= ScaledStep;
 TheParticle -> Position.Y -= ScaledStep;
 TheForce.Y = 0.5 * ( InteractionPotential() - TheForce.Y ) / ScaledStep; 
 TheParticle -> Position.Y += ScaledStep; 
 
  // calculate derivative along Z through threepoint finite difference... 
 ScaledStep = fabs( TheOldPosition.Z ) * DifferenceDelta;
 TheParticle -> Position.Z += ScaledStep;
 TheForce.Z = InteractionPotential();
 TheParticle -> Position.Z -= ScaledStep;
 TheParticle -> Position.Z -= ScaledStep;
 TheForce.Z = 0.5 * ( InteractionPotential() - TheForce.Z ) / ScaledStep; 
 TheParticle -> Position.Z += ScaledStep;
 
 // restore particles position...
 TheParticle -> Position = TheOldPosition;
   
 return( TheForce ); 
}

/* ----------------------------------------------------------------------- */

Vector FiniteDifferencePointPairInteraction::PointPairForce( Particle& P1, Particle& P2 )
{
 Vector TheForce, StepVector, OldPosition;
 double dV, dR;
 
 TheForce.X = TheForce.Y = TheForce.Z = 0.0;
 
 // check if we got the same thing twice...
 if( &(P1) == &(P2) ) return( TheForce );
  
 // get direction of the force on P1...
 TheForce = P1.Position - P2.Position;
 
 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
  TheForce = MimimumSystemImageConventionVector( TheForce, (*ThePeriodicBoxPointer) );
 
 dR = Abs( TheForce );
 TheForce = ( dR <= 0.0 ? TheForce.X = 1.0, TheForce : TheForce / dR );
 
 // get step size for finite difference...
 dR = DifferenceDelta * Abs( P1.Position );
 StepVector = TheForce * dR;
  
 // now do finite difference on particle P1 and calculate potential change...
 OldPosition = P1.Position;
 P1.Position += StepVector;
 dV = PointPairPotential( P1, P2 );
 P1.Position -= StepVector;
 P1.Position -= StepVector;
 dV = PointPairPotential( P1, P2 ) - dV;
 TheForce = TheForce * 0.5 * ( dV / dR );
 P1.Position = OldPosition;
 
 return( TheForce );
}

/* ----------------------------------------------------------------------- */

double FiniteDifferenceBendingInteraction::DerivativeOfBendingPotential( double cosangle )
{
 return( 0.5 * ( BendingPotential( cosangle + DifferenceDelta ) - BendingPotential( cosangle - DifferenceDelta ) ) / DifferenceDelta );
}

/* ----------------------------------------------------------------------- */

double FiniteDifferenceTorsionalInteraction::DerivativeOfTorsionalPotential( double cosangle )
{
  return( 0.5 * ( TorsionalPotential( cosangle + DifferenceDelta ) - TorsionalPotential( cosangle - DifferenceDelta ) ) / DifferenceDelta );
}

/* ----------------------------------------------------------------------- */