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

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

double HarmonicInteraction::PointPairPotential( Particle& P1, Particle& P2 )
{
 Vector dR;
 double r;
 
 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 ); 
 
 return( 0.5 * InteractionFactor * (r-R) * (r-R) ); 
}

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

Vector HarmonicInteraction::PointPairForce( Particle& P1, Particle& P2 )
{
 Vector TheForce;
 double r;
 
 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 );
 
 TheForce *= InteractionFactor * (r-R);
 return( TheForce );
}

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