/*   []----------------------------------------[]
      |               Torsion.cpp              |
     []----------------------------------------[]
      |                                        |
      |   AUTHOR:    MFSomers 2005.            |
      |   USE:       Bond torsion 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 "Torsional.h"

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

double TorsionalInteraction::InteractionPotential( void )
{
 int             Pos1, Pos2, Pos3, Pos4;
 double          TotalTorsionalPotential;
  
 TotalTorsionalPotential = 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;

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

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

Vector TorsionalInteraction::InteractionForce( unsigned int ID )
{
 int             Pos1, Pos2, Pos3, Pos4;
 Vector          TotalTorsionalVector;
  
 TotalTorsionalVector.X = 0.0;
 TotalTorsionalVector.Y = 0.0; 
 TotalTorsionalVector.Z = 0.0; 
 
 Pos1 = IsParticleIncluded( ID );
 
 if( Pos1 <= -1 ) RETURN_ERROR( "Particle not found", TotalTorsionalVector ); 
  
 // In general four bonds can contribute to the torsional force on a given particle.
 // 
 // Given the chain A-B-C-ID-D-E-F, the torsional angle A-B-C-ID (bond B-C), the angle
 // B-C-ID-D (bond C-ID), the angle C-ID-D-E (bond ID-D) and finally the torsional angle
 // ID-D-E-F (bond D-E) contribute to the force on particle ID.
 
 // first deal with the ID-D-E-F situation... find D E and F particles in chain ...
 Pos2 = Pos1 + 1;
 while( ( Pos2 < Particles.size() ) && ( Particles[ Pos2 ] == NULL ) ) ++Pos2;
 
 Pos3 = Pos2 + 1; 
 while( ( Pos3 < Particles.size() ) && ( Particles[ Pos3 ] == NULL ) ) ++Pos3;
 
 Pos4 = Pos3 + 1;
 while( ( Pos4 < Particles.size() ) && ( Particles[ Pos4 ] == NULL ) ) ++Pos4;
 
 if( Pos4 < Particles.size() )
  TotalTorsionalVector += TorsionalForce( *(Particles[ Pos1 ]), *(Particles[ Pos2 ]), *(Particles[ Pos3 ]), *(Particles[ Pos4 ]), 1 );
  
 // deal with the C-ID-D-E situation.. find particle C...
 Pos4 = Pos1 - 1;
 while( ( Pos4 >= 0 ) && ( Particles[ Pos4 ] == NULL ) ) --Pos4;
 
 if( ( Pos3 < Particles.size() ) && ( Pos4 >= 0 ) )
  TotalTorsionalVector += TorsionalForce( *(Particles[ Pos4 ]), *(Particles[ Pos1 ]), *(Particles[ Pos2 ]), *(Particles[ Pos3 ]), 2 );
 
 // deal with the B-C-ID-D situation... find particle B...
 Pos3 = Pos4 - 1;  
 while( ( Pos3 >= 0 ) && ( Particles[ Pos3 ] == NULL ) ) --Pos3;
  
 if( ( Pos2 < Particles.size() ) && ( Pos3 >= 0 ) )
  TotalTorsionalVector += TorsionalForce( *(Particles[ Pos3 ]), *(Particles[ Pos4 ]), *(Particles[ Pos1 ]), *(Particles[ Pos2 ]), 3 );
  
 // and finally deal with the A-B-C-ID situation... find particle A... 
 Pos2 = Pos3 - 1;
 while( ( Pos2 >= 0 ) && ( Particles[ Pos2 ] == NULL ) ) --Pos2;
    
 if( Pos2 >= 0 )
  TotalTorsionalVector += TorsionalForce( *(Particles[ Pos2 ]), *(Particles[ Pos3 ]), *(Particles[ Pos4 ]), *(Particles[ Pos1 ]), 4 ); 
 
 return( TotalTorsionalVector );
}

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

double TorsionalInteraction::TorsionalPotential( unsigned int ID1, unsigned int ID2, unsigned ID3, unsigned ID4 )
{
 int Pos1, Pos2, Pos3, Pos4;
 
 if( ( ID1 == ID2 ) || ( ID1 == ID3 ) || ( ID1 == ID4 ) || ( ID2 == ID3 ) || ( ID2 == ID4 ) || ( ID3 == ID4 ) ) 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 );  

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

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

Vector TorsionalInteraction::TorsionalForce( unsigned int ID1, unsigned int ID2, unsigned ID3, unsigned ID4, int p )
{
 int    Pos1, Pos2, Pos3, Pos4;
 Vector ZeroVector;
 
 ZeroVector.X = ZeroVector.Y = ZeroVector.Z = 0.0;
 
 if( ( p < 1 ) || ( p > 4 ) ) RETURN_ERROR( "Invalid p value", ZeroVector );
 
 if( ( ID1 == ID2 ) || ( ID1 == ID3 ) || ( ID1 == ID4 ) || ( ID2 == ID3 ) || ( ID2 == ID4 ) || ( ID3 == ID4 ) ) RETURN_ERROR( "7", 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 );  

 Pos4 = IsParticleIncluded( ID4 );
 if( Pos4 <= -1 ) RETURN_ERROR( "Particle not found", ZeroVector );    
    
 return( TorsionalForce( *(Particles[ Pos1 ]), *(Particles[ Pos2 ]), *(Particles[ Pos3 ]), *(Particles[ Pos4 ]), p ) );  
}

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

double TorsionalInteraction::TorsionalPotential( Particle& P1, Particle& P2, Particle& P3, Particle& P4 )
{
 Vector Rij, Rkj, Rkl, Rim, Rln;
 double cosangle, absRim, absRln, absRkjSQ;
 
 // calculate cosine of torsional angle first
 
 Rkj = P3.Position - P2.Position;
 
 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
  Rkj = MimimumSystemImageConventionVector( Rkj, (*ThePeriodicBoxPointer) );  
  
 absRkjSQ = Dot( Rkj, Rkj );
 if( absRkjSQ < CONST_EPSILON ) RETURN_ERROR( "Particles too close", TorsionalPotential( 0.0 ) );
 absRkjSQ = 1.0 / absRkjSQ;
  
 Rij = P1.Position - P2.Position;
 
 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
  Rij = MimimumSystemImageConventionVector( Rij, (*ThePeriodicBoxPointer) );  
  
 Rim = Rij - Dot( Rij, Rkj ) * Rkj * absRkjSQ;
 
 absRim = Abs( Rim );
 if( absRim < CONST_EPSILON ) return( TorsionalPotential( 0.0 ) );
 
 Rkl = P3.Position - P4.Position;
 
 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
  Rkl = MimimumSystemImageConventionVector( Rkl, (*ThePeriodicBoxPointer) );   
 
 Rln = Dot( Rkl, Rkj ) * Rkj * absRkjSQ - Rkl;
  
 absRln = Abs( Rln );
 if( absRln < CONST_EPSILON ) return( TorsionalPotential( 0.0 ) );
 
 absRim = 1.0 / absRim;
 absRln = 1.0 / absRln;
 
 cosangle = absRln * absRim * Dot( Rim, Rln );
 
 return( TorsionalPotential( cosangle ) );
}

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

Vector TorsionalInteraction::TorsionalForce( Particle& P1, Particle& P2, Particle& P3, Particle& P4, int p )
{
 Vector Rij, Rkj, Rkl, Rim, Rln, ForceVector;
 double cosangle, absRim, absRln, absRkjSQ, dVdcosangle;
 
 ForceVector.X = ForceVector.Y = ForceVector.Z = 0.0;
 
 if( ( p < 1 ) || ( p > 4 ) ) RETURN_ERROR( "Invalid p value", ForceVector );
 
 // calculate cosine of torsional angle first
 
 Rkj = P3.Position - P2.Position;

 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
  Rkj = MimimumSystemImageConventionVector( Rkj, (*ThePeriodicBoxPointer) );   
  
 absRkjSQ = Dot( Rkj, Rkj );
 if( absRkjSQ < CONST_EPSILON ) RETURN_ERROR( "Particles too close", ForceVector );
 absRkjSQ = 1.0 / absRkjSQ;
  
 Rij = P1.Position - P2.Position;

 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
  Rij = MimimumSystemImageConventionVector( Rij, (*ThePeriodicBoxPointer) );   
  
 Rim = Rij - Dot( Rij, Rkj ) * Rkj * absRkjSQ;
 
 absRim = Abs( Rim );
 if( absRim < CONST_EPSILON ) return( ForceVector );
 
 Rkl = P3.Position - P4.Position;
 
 // check if we have periodicity
 if( ThePeriodicBoxPointer != NULL )
  Rkl = MimimumSystemImageConventionVector( Rkl, (*ThePeriodicBoxPointer) );   
 
 Rln = Dot( Rkl, Rkj ) * Rkj * absRkjSQ - Rkl;
  
 absRln = Abs( Rln );
 if( absRln < CONST_EPSILON ) return( ForceVector );
 
 absRim = 1.0 / absRim;
 absRln = 1.0 / absRln;

 Rln = Rln * absRln;
 Rim = Rim * absRim;
  
 cosangle = Dot( Rim, Rln );
 dVdcosangle = DerivativeOfTorsionalPotential( cosangle );
 
 // we can have four different cases
 switch( p )
 {
  case 1: ForceVector = ( Rim * cosangle - Rln ) * absRim;
          break;
  case 2: ForceVector = ( Rim * cosangle - Rln ) * absRim * ( absRkjSQ * Dot( Rij, Rkj ) - 1 );
          ForceVector += ( Rim - Rln * cosangle ) * absRln * ( absRkjSQ * Dot( Rkl, Rkj ) ); 
          break;
  case 3: ForceVector = ( Rln - Rim * cosangle ) * absRim;
          ForceVector += ( Rim - Rln * cosangle ) * absRln;
	  ForceVector += ( Rln - Rim * cosangle ) * absRim * ( absRkjSQ * Dot( Rij, Rkj ) - 1 );
	  ForceVector += ( Rln  * cosangle - Rim ) * absRln * ( absRkjSQ * Dot( Rkl, Rkj ) ); 
          break;
  case 4: ForceVector = ( Rln * cosangle - Rim ) * absRln;
          break;
 }
 
 return( dVdcosangle * ForceVector );
}

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

double TorsionalInteraction::TorsionalPotential( double cosangle )
{
 return( 0.0 );
}

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

double TorsionalInteraction::DerivativeOfTorsionalPotential( double cosangle )
{
 return( 0.0 );
}

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

double PeriodicTorsionalInteraction::TorsionalPotential( double cosangle )
{
 if( fabs( cosangle ) > 1.0 - CONST_EPSILON ) 
  return( InteractionFactor * ( 1 - cos( N * Angle ) ) );
 else
  return( InteractionFactor * ( 1 - cos( N * ( acos( cosangle ) - Angle ) ) ) );
}

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

double PeriodicTorsionalInteraction::DerivativeOfTorsionalPotential( double cosangle )
{
 double X, Y;
 
 if( fabs( cosangle ) > 1.0 - CONST_EPSILON ) 
  X = InteractionFactor * N * sin( N * Angle );
 else
  X = InteractionFactor * N * sin( N * ( Angle - acos( cosangle ) ) );

 Y = 1.0 - cosangle * cosangle;

 if( Y < CONST_EPSILON ) return( 0.0 );
  
 return( X / sqrt( Y ) );
}

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