/*   []----------------------------------------[]
      |              Bending.cpp               |
     []----------------------------------------[]
      |                                        |
      |   AUTHOR:    MFSomers 2005.            |
      |   USE:       Bond bending interaction  |
      |               functions...             |      
      |                                        |
     []----------------------------------------[]
*/

// 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 "Bending.h"

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

double BendingInteraction::InteractionPotential( void )
{
 int             Pos1, Pos2, Pos3;
 double          TotalBendingPotential;
 
 TotalBendingPotential = 0.0;
 
 // find first particle in list
 Pos1 = 0;
 while( ( Pos1 < Particles.size() ) && ( Particles[ Pos1 ] == NULL ) ) ++Pos1;
 
 // find second particle 
 Pos2 = Pos1 + 1;
 while( ( Pos2 < Particles.size() ) && ( Particles[ Pos2 ] == NULL ) ) ++Pos2;

 // and find third particle 
 Pos3 = Pos2 + 1;
 while( ( Pos3 < Particles.size() ) && ( Particles[ Pos3 ] == NULL ) ) ++Pos3;
 
 // now loop along chain and add potential contribution of all angles...
 while( Pos3 < Particles.size() )
 {
  TotalBendingPotential += BendingPotential( *(Particles[ Pos1 ]), *(Particles[ Pos2 ]), *(Particles[ Pos3 ]) ); 
 
  Pos1 = Pos2;
  Pos2 = Pos3;
  ++Pos3;
  while( ( Pos3 < Particles.size() ) && ( Particles[ Pos3 ] == NULL ) ) ++Pos3;
 }
 
 return( TotalBendingPotential ); 
} 

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

Vector BendingInteraction::InteractionForce( unsigned int ID )
{
 int             Pos1, Pos2, Pos3;
 Vector          TotalBendingVector;
 
 TotalBendingVector.X = 0.0;
 TotalBendingVector.Y = 0.0; 
 TotalBendingVector.Z = 0.0; 
 
 Pos1 = IsParticleIncluded( ID );
 
 if( Pos1 <= -1 ) RETURN_ERROR( "Particle not found", TotalBendingVector ); 
  
 // at most three angles in the chain can contribute to the force on the particle of given ID.
 //
 // given chain A-B-ID-C-D:
 //
 // angle ID-C-D, angle B-ID-C and angle A-B-ID all contribute to the force on particle ID
 
 // find second particle further in chain 
 Pos2 = Pos1 + 1;
 while( ( Pos2 < Particles.size() ) && ( Particles[ Pos2 ] == NULL ) ) ++Pos2;
 
 // find third particle further in chain
 Pos3 = Pos2 + 1;
 while( ( Pos3 < Particles.size() ) && ( Particles[ Pos3 ] == NULL ) ) ++Pos3;
 
 // add to vector if angle ID-C-D is defined
 if( Pos3 < Particles.size() )
  TotalBendingVector += BendingForce( *(Particles[ Pos1 ]), *(Particles[ Pos2 ]), *(Particles[ Pos3 ]), 1 );
  
 // now find first particle back in chain
 Pos3 = Pos1 - 1;
 while( ( Pos3 >= 0 ) && ( Particles[ Pos3 ] == NULL ) ) --Pos3;  
 
 // if angle B-ID-C is present, add to vector
 if( ( Pos3 >= 0 ) && ( Pos2 < Particles.size() ) )
  TotalBendingVector += BendingForce( *(Particles[ Pos3 ]), *(Particles[ Pos1 ]), *(Particles[ Pos2 ]), 2 );
  
 // finally find second particle back in chain
 Pos2 = Pos3 - 1;
 while( ( Pos2 >= 0 ) && ( Particles[ Pos2 ] == NULL ) ) --Pos2;
   
 // if angle A-B-ID is defined, add to vector 
 if( Pos2 >= 0 ) 
  TotalBendingVector += BendingForce( *(Particles[ Pos2 ]), *(Particles[ Pos3 ]), *(Particles[ Pos1 ]), 3 );

 return( TotalBendingVector );
}

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

double BendingInteraction::BendingPotential( unsigned int ID1, unsigned int ID2, unsigned ID3 )
{
 int Pos1, Pos2, Pos3;
 
 if( ( ID1 == ID2 ) || ( ID2 == ID3 ) || ( ID3 == ID1 ) ) RETURN_ERROR( "Two or more identical particles", 0.0 );
 
 Pos1 = IsParticleIncluded( ID1 );
 if( Pos1 <= -1 ) RETURN_ERROR( "Particle not found", 0.0 );
  
 Pos2 = IsParticleIncluded( ID2 );
 if( Pos2 <= -1 ) RETURN_ERROR( "Particle not found", 0.0 );

 Pos3 = IsParticleIncluded( ID3 );
 if( Pos3 <= -1 ) RETURN_ERROR( "Particle not found", 0.0 );  
    
 return( BendingPotential( *(Particles[ Pos1 ]), *(Particles[ Pos2 ]), *(Particles[ Pos3 ]) ) );  
}

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

Vector BendingInteraction::BendingForce( unsigned int ID1, unsigned int ID2, unsigned ID3, int p )
{
 int    Pos1, Pos2, Pos3;
 Vector ZeroVector;
 
 ZeroVector.X = ZeroVector.Y = ZeroVector.Z = 0.0;
 
 if( ( p > 3 ) || ( p < 1 ) ) RETURN_ERROR( "Invalid p value", ZeroVector );
 
 if( ( ID1 == ID2 ) || ( ID2 == ID3 ) || ( ID3 == ID1 ) ) RETURN_ERROR( "Two or more identical particles", ZeroVector );
 
 Pos1 = IsParticleIncluded( ID1 );
 if( Pos1 <= -1 ) RETURN_ERROR( "Particle not found", ZeroVector );
  
 Pos2 = IsParticleIncluded( ID2 );
 if( Pos2 <= -1 ) RETURN_ERROR( "Particle not found", ZeroVector );
  
 Pos3 = IsParticleIncluded( ID3 );
 if( Pos3 <= -1 ) RETURN_ERROR( "Particle not found", ZeroVector );  
  
 return( BendingForce( *(Particles[ Pos1 ]), *(Particles[ Pos2 ]), *(Particles[ Pos3 ]), p ) );  
}

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

double BendingInteraction::BendingPotential( Particle& P1, Particle& P2, Particle& P3 )
{
 Vector Rij, Rkj;
 double cosangle, absRij, absRkj;
 
 cosangle = 1.0;
 
 Rij = P1.Position - P2.Position;
 Rkj = P3.Position - P2.Position;
 
 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
 {
  Rij = MimimumSystemImageConventionVector( Rij, (*ThePeriodicBoxPointer) );
  Rkj = MimimumSystemImageConventionVector( Rkj, (*ThePeriodicBoxPointer) );
 }
 
 absRij = Abs( Rij );
 absRkj = Abs( Rkj );

 if( absRij < CONST_EPSILON ) RETURN_ERROR( "Particles too close", 0.0 );
 if( absRkj < CONST_EPSILON ) RETURN_ERROR( "Particles too close", 0.0 );
 
 if( ( absRkj >= CONST_EPSILON ) && ( absRij >= CONST_EPSILON ) )
  cosangle = Dot( Rij, Rkj ) / ( absRij * absRkj );
   
 return( BendingPotential( cosangle ) );
}

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

Vector BendingInteraction::BendingForce( Particle& P1, Particle& P2, Particle& P3, int p )
{
 Vector ForceVector;
 Vector Rij, Rkj;
 double cosangle, absRij, absRkj, dVdcosangle;
 
 ForceVector.X = ForceVector.Y = ForceVector.Z = 0.0;
 
 if( ( p > 3 ) || ( p < 1 ) ) RETURN_ERROR( "Invalid p value", ForceVector );
 
 Rij = P1.Position - P2.Position;
 Rkj = P3.Position - P2.Position;
 
 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
 {
  Rij = MimimumSystemImageConventionVector( Rij, (*ThePeriodicBoxPointer) );
  Rkj = MimimumSystemImageConventionVector( Rkj, (*ThePeriodicBoxPointer) );
 } 
 
 absRij = Abs( Rij );
 if( absRij < CONST_EPSILON ) RETURN_ERROR( "Particles too close", ForceVector );
 absRij = 1.0 / absRij;
 
 absRkj = Abs( Rkj );
 if( absRkj < CONST_EPSILON ) RETURN_ERROR( "Particles too close", ForceVector );
 absRkj = 1.0 / absRkj;
 
 Rij = absRij * Rij;
 Rkj = absRkj * Rkj;
 
 cosangle = Dot( Rij, Rkj );
 dVdcosangle = DerivativeOfBendingPotential( cosangle );
 
 // we can have three cases, the force on particle 1, 2 or 3 is requested
 
 switch( p )
 {
  case 1: ForceVector = ( Rij * cosangle - Rkj ) * absRij; 
          break;
  case 2: ForceVector = ( Rkj - Rij * cosangle ) * absRij;
          ForceVector += ( Rij - Rkj * cosangle ) * absRkj;
          break;
  case 3: ForceVector = ( Rkj * cosangle - Rij ) * absRkj;
          break;
 }
 
 return( dVdcosangle * ForceVector );
}

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

double BendingInteraction::BendingPotential( double cosangle )
{
 return( 0.0 );
}

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

double BendingInteraction::DerivativeOfBendingPotential( double cosangle )
{
 return( 0.0 );
}

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

double HarmonicBendingInteraction::BendingPotential( double cosangle )
{
 return( 0.5 * InteractionFactor * ( cosangle - CosAngle ) * ( cosangle - CosAngle ) );
}

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

double HarmonicBendingInteraction::DerivativeOfBendingPotential( double cosangle )
{
 return( InteractionFactor * ( cosangle - CosAngle ) );
}

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