/*   []----------------------------------------[]
      |            LennardJones.cpp            |
     []----------------------------------------[]
      |                                        |
      |   AUTHOR:    MFSomers 2005.            |
      |   USE:       Defines a LennardJones    |
      |              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 "LennardJones.h"

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

double LennardJonesInteraction::PointPairPotential( Particle& P1, Particle& P2 )
{
 Vector dR;
 double r, a;
 
 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 ); 
 
 if( r < CONST_EPSILON ) RETURN_ERROR( "Particles too close", CONST_HUGE_VALUE );
  
 a = ( R / r );
 a = a * a;
 a = a * a * a;             /* a = (R/r) ^ 6 */
 
 return( Epsilon * a * ( a - 2.0 ) ); 
}

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

Vector LennardJonesInteraction::PointPairForce( Particle& P1, Particle& P2 )
{
 Vector TheForce;
 double r, a;
 
 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 = CONST_HUGE_VALUE, TheForce : TheForce / r );
 if( r < CONST_EPSILON ) RETURN_ERROR( "Particles too close", TheForce );
  
 a = ( R / r );
 a = a * a;
 a = a * a * a;             /* a = (R/r) ^ 6 */
 
 TheForce *= ( 12.0 * Epsilon * a * ( 1 - a ) / r );
 
 return( TheForce );
}

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