/*   []----------------------------------------[]
      |                Morse.cpp               |
     []----------------------------------------[]
      |                                        |
      |   AUTHOR:    MFSomers 2005.            |
      |   USE:       Defines a Morse type      |
      |               interaction...           |
      |                                        |
     []----------------------------------------[]
*/

// 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 "Morse.h"

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

double MorseInteraction::PointPairPotential( Particle& P1, Particle& P2 )
{
 Vector dR;
 double r, exp_part;
 
 if( &(P1) == &(P2) ) RETURN_ERROR( "Identical particles", 0.0 );
 
 dR = P1.Position - P2.Position;
 
 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
  dR = MimimumSystemImageConventionVector( dR, (*ThePeriodicBoxPointer) ); 
  
  // determine distance between particles ...  
 r = Abs( dR ); 
 
 exp_part = 1 - exp( Alfa * ( R - r ) );
  
 return( D * exp_part * exp_part - D ); 
}

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

Vector MorseInteraction::PointPairForce( Particle& P1, Particle& P2 )
{
 Vector TheForce;
 double r, exp_part;
 
 TheForce.X = TheForce.Y = TheForce.Z = 0.0;
 
 if( &(P1) == &(P2) ) RETURN_ERROR( "Identical particles", TheForce );
  
 // get distance between particles and direction of force... 
 TheForce = P2.Position - P1.Position;
 
 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
  TheForce = MimimumSystemImageConventionVector( TheForce, (*ThePeriodicBoxPointer) );  
 
 r = Abs( TheForce ); 
 TheForce = ( r < CONST_EPSILON ? TheForce.X = 1.0, TheForce : TheForce / r );
 
 exp_part = exp( Alfa * ( R - r ) );
 
 TheForce *= 2.0 * Alfa * D * exp_part * ( 1 - exp_part );
 
 return( TheForce );
}

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