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

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

double CoulombInteraction::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 ); 

 if( r < CONST_EPSILON ) RETURN_ERROR( "Particles too close", CONST_HUGE_VALUE );

 return( InteractionFactor * P1.Charge * P2.Charge / r ); 
}

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

Vector CoulombInteraction::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 = P1.Position - P2.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 ); 

 if( r < CONST_EPSILON ) RETURN_ERROR( "Particles too close", TheForce );

 TheForce = InteractionFactor * P1.Charge * P2.Charge * TheForce / (r*r*r);
 
 return( TheForce );
}

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