/*   []----------------------------------------[]
      |            SteepestDescent.h           |
     []----------------------------------------[]
      |                                        |
      |   AUTHOR:    MFSomers 2005.            |
      |   USE:       Implements a conformation |
      |              searcher using steepest   |
      |              descent method...         |
      |                                        |
     []----------------------------------------[]
*/

// 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 "SteepestDescent.h"
#include "DistanceConstraint.h"

/* -------------------------------------------------- SteepestDescent ---- */

int SteepestDescent::FindMinima( int MaxSteps, double MaxStepSize, double Accuracy, ConstraintListPointer Constraints )
{
 DistanceConstraintPointer C;
 double StepSize, Change;
 int n;
 
 if( MaxSteps <= 0 ) return( 0 );
  
 // first we setup the constraint tables indices to particles once at the start of a steepest descent search
 
 if( Constraints != NULL )
  for( int CIndex = 0; CIndex < Constraints -> size(); ++CIndex )
   if( ( C = (DistanceConstraintPointer)((*Constraints)[ CIndex ]) ) != NULL )
    if( ( C -> Type == AGeneralDistanceConstraint ) && ( C -> Distance > 0.0 ) )
     for( int n = 0; n < C -> GetParticleListSize(); ++n )   
      if( C -> ParticleIDs[ n ] ) 
       C -> ParticleIndices[ n ] = IsParticleIncluded( C -> ParticleIDs[ n ] );
 
 Accuracy = fabs( Accuracy );
 Accuracy = Accuracy >= CONST_EPSILON ? Accuracy : CONST_EPSILON;
 
 MaxStepSize = fabs( MaxStepSize );
 MaxStepSize = MaxStepSize >= CONST_EPSILON ? MaxStepSize : CONST_EPSILON;
 
 StepSize = MaxStepSize;

 for( n = 0; n < MaxSteps; ++n )                          // step a given max number of times
 {
  if( !StepSteepestDescent( StepSize, Constraints ) ) break;     // stop if positions didn't change

  Change = RelativeTotalInteractionPotentialChange();

  if( fabs( Change ) <= Accuracy )                        // now check if converging to accuracy
  {                                                       // both in the change of the potential 
   if( MaxStepSize <= Accuracy ) break;                      // and in the step size
   StepSize *= 0.8;
   MaxStepSize *= 0.8;
  }
  
  CycleTheStates( -1 );                                   // cycle the states for the next step  
    
  if( Change > 0.0 )                                      // and change step size accordingly
  {
   StepSize *= 0.5;
   if( StepSize < CONST_EPSILON ) StepSize = CONST_EPSILON;
  }
  else
  {
   StepSize *= 1.2;
   if( StepSize >= MaxStepSize ) StepSize = MaxStepSize;
  } 
  
 }
 
 return( n );
}

/* ----------------------------------------------------------------------- */
/* This algorithm combines the SHAKE algorithm with a step towards the     */
/* forces. See for more details on the SHAKE constraint algorith:          */
/*                                                                         */
/* H. C. Andersen, Journal of Computational Physics 52, 24-34 (1983)       */
/* http://www.CCL.net/cca/software/SOURCES/FORTRAN/allen-tildesley-book/f.08.shtml */

int SteepestDescent::StepSteepestDescent( double StepSize, ConstraintListPointer Constraints )
{
 ParticleListPointer CurrentState = TheParticleListList[ 1 ];
 ParticleListPointer NewState = TheParticleListList[ 0 ];
 int Done, ReturnFlag;
  
 if( CurrentState == NULL ) return ( 0 );
 if( NewState == NULL ) return( 0 );
 
 ReturnFlag = 0;

 #pragma omp parallel
 { 
  DistanceConstraintPointer C;                          // all local private variables
  ParticlePointer           TheParticle, TheSecondParticle;
  int                       n, TheParticleIndex, TheSecondParticleIndex, Size;
  Vector                    v, dv;    

  // first step towards direction of forces

  Size = CurrentState -> size();

  #pragma omp for schedule( static )
  for( n = 0; n < Size; ++n )
   if( ( ( TheParticle = (*CurrentState)[ n ] ) != NULL ) && ( ( TheSecondParticle = (*NewState)[ n ] ) != NULL ) )
   {
    v = TheParticle -> TheTotalForce();
    double Absv = Abs( v );
   
    if( Absv > 0.0 )
    {
     ReturnFlag = 1;
     TheSecondParticle -> Position = TheParticle -> Position + StepSize * v / Absv;
    
     // also mark the fact that the particle has been moved for constraint checks 
     TheSecondParticle -> TheInternalFlag = 1;    
    }
    else
     TheSecondParticle -> Position = TheParticle -> Position;
   }
 
  // now apply possible distance constraints through SHAKE
 
  if( Constraints != NULL )
  {
   int Itteration = 0;
   
   Size = Constraints -> size();

   do
   {
    #pragma omp barrier

    Done = 1; 

    #pragma omp for schedule( static )
    for( int CIndex = 0; CIndex < Size; ++CIndex )
     if( ( C = (DistanceConstraintPointer)((*Constraints)[ CIndex ]) ) != NULL )
      if( ( C -> Type == AGeneralDistanceConstraint ) && ( C -> Distance > 0.0 ) && ( Itteration < C -> MaxItterations ) )
      {
       double DistanceSqr = ( C -> Distance ) * ( C -> Distance );
          
       for( int PIndex = n = 0; n < C -> NrOfParticlesIncluded() - 1; ++n )
       {
        while( ( TheParticleIndex = C -> ParticleIndices[ PIndex ] ) < 0 ) ++PIndex;          // get first and second particle 
        if( (*CurrentState)[ TheParticleIndex ]-> Mass == 0.0 ) continue;                     // and make sure they have a mass
       
        ++PIndex; while( ( TheSecondParticleIndex = C -> ParticleIndices[ PIndex ] ) < 0 ) ++PIndex;
        if( (*CurrentState)[ TheSecondParticleIndex ]-> Mass == 0.0 ) continue;
       
        TheParticle       = (*NewState)[ TheParticleIndex ];
        TheSecondParticle = (*NewState)[ TheSecondParticleIndex ];
       
        // if one of these particles has been moved because of constraints or whatever
        if( ( TheParticle -> TheInternalFlag > 0 ) || ( TheSecondParticle -> TheInternalFlag > 0 ) )
        {
 
         Vector PAB = TheParticle -> Position - TheSecondParticle -> Position;
 	 double AbsPABSqr = AbsSqr( PAB );
	 double DiffSqr = DistanceSqr - AbsPABSqr;
	
	 if( fabs( DiffSqr ) > ( 2.0 * (C -> Tolerance) * DistanceSqr ) )
	 {
		
     	  Vector RAB = (*CurrentState)[ TheParticleIndex ] -> Position - (*CurrentState)[ TheSecondParticleIndex ] -> Position;
	  double RABDotPAB = Dot( RAB, PAB );
	 
	  double InvMassA = ( 1.0 / TheParticle -> Mass );
	  double InvMassB = ( 1.0 / TheSecondParticle -> Mass );
	 
	  double GAB = 0.5 * DiffSqr / ( ( InvMassA + InvMassB ) * RABDotPAB );
	 
	  dv = GAB * RAB;
	 
	  TheParticle -> Position       += dv * InvMassA;
	  TheSecondParticle -> Position -= dv * InvMassB;

	  // signal we need to continue	 
	  Done = 0;
	 
	  // and that these two particles have been moved so convergence checks are again performed
	  TheParticle -> TheInternalFlag = 2;
	  TheSecondParticle -> TheInternalFlag = 2;
	 
	  // and the steepest descent steps are not yet finished cause particles moved
	  ReturnFlag = 1;
	 } 
	
	
 	 TheParticle -> TheInternalFlag--;
	 TheSecondParticle -> TheInternalFlag--;
        }
       
       }
    
      }
   
    Itteration++;    

    #pragma omp barrier
   } while( !Done );
  
  } 
  
 } // end parallel region
 
 return( ReturnFlag );   
}

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

double SteepestDescent::RelativeTotalInteractionPotentialChange( void )
{
 InteractionListPointer CurrentState = TheInteractionListList[ 1 ];
 InteractionListPointer NewState = TheInteractionListList[ 0 ];
 double CurrentEnergy, NewEnergy;
 
 if( NewState != NULL && CurrentState != NULL )
 {
  CurrentEnergy = TotalPotentialEnergy( (*CurrentState) );
  NewEnergy     = TotalPotentialEnergy( (*NewState) );
  
  if( NewEnergy != 0.0 )
   return( ( NewEnergy - CurrentEnergy ) / fabs( NewEnergy ) );   
  else 
   return( NewEnergy - CurrentEnergy );
 } 
 else
  return( 0.0 ); 
}

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