// A simple LennardJones force field...

#if defined(__OMP__)
#include <omp.h>
#endif

#include "LennardJones.h"

// ------------------------------------------------------------

int LennardJonesForces( int N, Particle P[] )
{
 Vector r;
 int i, j;
 double k, kk;
 
#if defined(__OMP__)

#pragma omp parallel for private(i,j,r,k,kk)
 for( i = 0; i < N; ++i )
 {
  P[ i ].F.X = 0.0;
  P[ i ].F.Y = 0.0;
  P[ i ].F.Z = 0.0;

  for( j = 0; j < N; ++j )
   if( i != j )
   {
    r = SubVectors( P[ j ].r, P[ i ].r );
    k = 1.0 / AbsVector( r );
    kk = LENNARDJONES_SIGMA * k;
    kk = kk * kk * kk;
    kk = kk * kk;
    k = 4.0 * LENNARDJONES_EPSILON * k * ( 12.0 * kk * kk - 6.0 * kk );
    P[ i ].F = AddVectors( P[ i ].F, ScalarVector( k, r ) );
   }
 }

#else

 for( i = 0; i < N; ++i )
 {
  P[ i ].F.X = 0.0;
  P[ i ].F.Y = 0.0;
  P[ i ].F.Z = 0.0;
 }

 for( i = 0; i < N; ++i )
  for( j = i + 1; j < N; ++j )
  {
   r = SubVectors( P[ j ].r, P[ i ].r );
#  if defined(__FAST__)
#  warning Using profiled code...
   k = 1.0 / AbsVector( r );
   kk = LENNARDJONES_SIGMA * k;
   kk = kk * kk * kk;
   kk = kk * kk;
   k = 4.0 * LENNARDJONES_EPSILON * k * ( 12.0 * kk * kk - 6.0 * kk );
#  else
   k = 4.0 * LENNARDJONES_EPSILON * ( 12.0 * pow( LENNARDJONES_SIGMA / AbsVector( r ), 12.0 ) - 6.0 * pow( LENNARDJONES_SIGMA / AbsVector( r ), 6.0 ) ) / AbsVector( r );
#  endif
   P[ i ].F = AddVectors( P[ i ].F, ScalarVector( k, r ) );
   P[ j ].F = SubVectors( P[ j ].F, ScalarVector( k, r ) );
  }

#endif

 return( 0 );
}

// ------------------------------------------------------------

double LennardJonesPotential( int N, Particle P[] )
{
 double r, V = 0.0;
 int i, j;
 
 for( i = 0; i < N; ++i )
  for( j = i + 1; j < N; ++j )
  {
   r = LENNARDJONES_SIGMA / AbsVector( SubVectors( P[ j ].r, P[ i ].r ) );
   V += 4.0 * LENNARDJONES_EPSILON * ( pow( r, 12.0 ) - pow( r, 6.0 ) );
  }

 return( V );
}

// ------------------------------------------------------------