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

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

double RydbergInteraction::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 ); 
 
 r  = r - R;  

 return( D - D * ( 1.0 + A1*r + A2*r*r + A3*r*r*r ) * exp( -1.0 * Gamma * r ) ); 
}

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

Vector RydbergInteraction::PointPairForce( Particle& P1, Particle& P2 )
{
 Vector TheForce;
 double r, r2, 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 );
 
 r  = r - R;
 r2 = r * r;
 exp_part = exp( -1.0 * Gamma * r );
 
 TheForce *= D * exp_part * ( (Gamma-A1) + (Gamma*A1-2.0*A2)*r + (Gamma*A2-3.0*A3)*r2 + Gamma*A3*r2*r );
 
 return( TheForce );
}

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