/*   []----------------------------------------[]
      |               ReadInput.cpp            |
     []----------------------------------------[]
      |                                        |
      |   AUTHOR:    MFSomers 2005.            |
      |   USE:       Read input files ...      |
      |                                        |
     []----------------------------------------[]
*/

// 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 "Coulomb.h"
#include "Gravity.h"
#include "Harmonic.h"
#include "Morse.h"
#include "LennardJones.h"
#include "Rydberg.h"

#include "NewtonDynamics.h"
#include "SteepestDescent.h"

#include "Bending.h"
#include "Torsional.h"

#include "DistanceConstraint.h"

#include "ReadInput.h"

#include <cstdio>
#include <cctype>
#include <cstring>

char __ATempStringForSPrintfs[ 1024 ];

/* ------------------------------ The interal string handler routines ---------- */

char *SearchAndSkipField( char *Line, char *Field )
{
 if( Line == NULL ) return( NULL );
 if( Field == NULL ) return( NULL );
 if( (*Field) == '\0' ) return( NULL );
 
 while( (*Line) != '\0' )
 {
  while( isspace( ( int )(*Line) ) ) ++Line;         // ignore spaces at start
  
  if( !strncmp( Line, Field, strlen( Field ) ) )     // found the field string
  {
   Line += strlen( Field );                          // and go beyond it to the rest of the line or '\0' probably
   while( isspace( ( int )(*Line) ) ) ++Line;        // then skip spaces and try to find the = char
   if( (*Line) == '=' ) return( Line + 1 );      
  }
  else
   while( !isspace( ( int )(*Line) ) && ( (*Line) != '\0' ) ) ++Line;        // run through the non-spaces or stop at EOS
   
 } 
 
 return( NULL );
}

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

int FieldPresent( char *Line, char *Field )
{
 return( SearchAndSkipField( Line, Field ) != NULL );
}

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

int GetInteger( char *Line, char *Field, int& I )
{
 if( ( Line = SearchAndSkipField( Line, Field ) ) == NULL ) return( 0 );
 return( sscanf( Line, " %d ", &I ) == 1 );
}

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

int GetInteger( char *Line, char *Field, unsigned int& I )
{
 if( ( Line = SearchAndSkipField( Line, Field ) ) == NULL ) return( 0 );
 while( isspace( ( int )(*Line) ) ) ++Line;
 if( !isdigit( ( int )(*Line) ) )  return( 0 );
 return( sscanf( Line, " %u ", &I ) == 1 );
}

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

int GetDouble( char *Line, char *Field, double& D )
{
 if( ( Line = SearchAndSkipField( Line, Field ) ) == NULL ) return( 0 );
 return( sscanf( Line, " %lf ", &D ) == 1 );
}

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

int GetVector( char *Line, char *Field, Vector& V )
{
 if( ( Line = SearchAndSkipField( Line, Field ) ) == NULL ) return( 0 );
 return( sscanf( Line, " [ %lf , %lf , %lf ] ", &V.X, &V.Y, &V.Z ) == 3 );
}

/* -------------------------------------------- The internal table routines ----- */ 

InteractionPointer ReadTheCoulombInteraction( char *LinePos, InputFileDataType& TheData )
{
 double factor = COULOMB_VACUUM_CONSTANT;
 
 if( FieldPresent( LinePos, "f" ) )
  if( !GetDouble( LinePos, "f", factor ) ) return( NULL );
  
 if( factor == 0.0 ) return( NULL );
 
 return( new CoulombInteraction( factor, TheData.DefaultParticle.Id, TheData.TheBoxPointer ) );
} 

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

char *StringTheCoulombInteraction( void *TheInteraction )
{
 CoulombInteractionPointer TheI = ( CoulombInteractionPointer )( TheInteraction );
 double f = TheI -> InteractionFactor;
 
 sprintf( __ATempStringForSPrintfs, " f=%.20le ", f );
 return( __ATempStringForSPrintfs );
}

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

InteractionPointer ReadTheGravitationalInteraction( char *LinePos, InputFileDataType& TheData )
{
 double factor = CONST_GRAVITATIONAL_CONST;
 
 if( FieldPresent( LinePos, "f" ) )
  if( !GetDouble( LinePos, "f", factor ) ) return( NULL );
 
 if( factor == 0.0 ) return( NULL ); 
  
 return( new GravitationalInteraction( factor, TheData.DefaultParticle.Id, TheData.TheBoxPointer ) );
} 

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

char *StringTheGravitationalInteraction( void *TheInteraction )
{
 GravitationalInteractionPointer TheI = ( GravitationalInteractionPointer )( TheInteraction );
 double f = TheI -> InteractionFactor;
 
 sprintf( __ATempStringForSPrintfs, " f=%.20le ", f );
 return( __ATempStringForSPrintfs );
}

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

InteractionPointer ReadTheHarmonicInteraction( char *LinePos, InputFileDataType& TheData )
{
 double factor, equilibrium;
 
 if( !GetDouble( LinePos, "f", factor ) ) return( NULL );
 if( !GetDouble( LinePos, "r0", equilibrium ) ) return( NULL ); 
   
 if( factor == 0.0 ) return( NULL ); 
  
 return( new HarmonicInteraction( factor, equilibrium, TheData.DefaultParticle.Id, TheData.TheBoxPointer ) );
} 

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

char *StringTheHarmonicInteraction( void *TheInteraction )
{
 HarmonicInteractionPointer TheI = ( HarmonicInteractionPointer )( TheInteraction );
 double f = TheI -> InteractionFactor;
 double r0 = TheI -> R;
 
 sprintf( __ATempStringForSPrintfs, " f=%.20le r0=%.20le ", f, r0 );
 return( __ATempStringForSPrintfs );
}

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

InteractionPointer ReadTheRydbergInteraction( char *LinePos, InputFileDataType& TheData )
{
 double factor, equilibrium, exponent;
 double a1 = 0.0;
 double a2 = 0.0;
 double a3 = 0.0;
 
 if( !GetDouble( LinePos, "f", factor ) ) return( NULL );
 if( !GetDouble( LinePos, "r0", equilibrium ) ) return( NULL ); 
 if( !GetDouble( LinePos, "exp", exponent ) ) return( NULL );   

 if( FieldPresent( LinePos, "a1" ) )
  if( !GetDouble( LinePos, "a1", a1 ) ) return( NULL );     
    
 if( FieldPresent( LinePos, "a2" ) )
  if( !GetDouble( LinePos, "a2", a2 ) ) return( NULL );     

 if( FieldPresent( LinePos, "a3" ) )
  if( !GetDouble( LinePos, "a3", a3 ) ) return( NULL );        
  
 if( factor == 0.0 ) return( NULL ); 
  
 return( new RydbergInteraction( factor, equilibrium, a1, a2, a3, exponent, TheData.DefaultParticle.Id, TheData.TheBoxPointer ) );
} 

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

char *StringTheRydbergInteraction( void *TheInteraction )
{
 RydbergInteractionPointer TheI = ( RydbergInteractionPointer )( TheInteraction );
 double f = TheI -> D;
 double r0 = TheI -> R;
 double exp = TheI -> Gamma;
 double a1 = TheI -> A1;
 double a2 = TheI -> A2;
 double a3 = TheI -> A3;
 
 sprintf( __ATempStringForSPrintfs, " f=%.20le r0=%.20le exp=%.20le a1=%.20le a2=%.20le a3=%.20le ", f, r0, exp, a1, a2, a3 );
 return( __ATempStringForSPrintfs );
}

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

InteractionPointer ReadTheMorseInteraction( char *LinePos, InputFileDataType& TheData )
{
 double factor, equilibrium, exponent;
 
 if( !GetDouble( LinePos, "f", factor ) ) return( NULL );
 if( !GetDouble( LinePos, "r0", equilibrium ) ) return( NULL ); 
 if( !GetDouble( LinePos, "exp", exponent ) ) return( NULL );   
   
 if( factor == 0.0 ) return( NULL ); 
 if( exponent == 0.0 ) return( NULL ); 
  
 return( new MorseInteraction( factor, equilibrium, exponent, TheData.DefaultParticle.Id, TheData.TheBoxPointer ) );
}

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

char *StringTheMorseInteraction( void *TheInteraction )
{
 MorseInteractionPointer TheI = ( MorseInteractionPointer )( TheInteraction );
 double f = TheI -> D;
 double r0 = TheI -> R;
 double exp = TheI -> Alfa;
 
 sprintf( __ATempStringForSPrintfs, " f=%.20le r0=%.20le exp=%.20le ", f, r0, exp );
 return( __ATempStringForSPrintfs );
}

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

InteractionPointer ReadTheLennardJonesInteraction( char *LinePos, InputFileDataType& TheData )
{
 double factor, equilibrium;
 
 if( !GetDouble( LinePos, "f", factor ) ) return( NULL );
 if( !GetDouble( LinePos, "r0", equilibrium ) ) return( NULL ); 
  
 if( factor == 0.0 ) return( NULL ); 
  
 return( new LennardJonesInteraction( factor, equilibrium, TheData.DefaultParticle.Id, TheData.TheBoxPointer ) );
}

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

char *StringTheLennardJonesInteraction( void *TheInteraction )
{
 LennardJonesInteractionPointer TheI = ( LennardJonesInteractionPointer )( TheInteraction );
 double f = TheI -> Epsilon;
 double r0 = TheI -> R;
  
 sprintf( __ATempStringForSPrintfs, " f=%.20le r0=%.20le ", f, r0 );
 return( __ATempStringForSPrintfs );
}

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

InteractionPointer ReadTheHarmonicBendingInteraction( char *LinePos, InputFileDataType& TheData )
{
 double factor, angle;
 
 if( !FieldPresent( LinePos, "rad" ) && !FieldPresent( LinePos, "deg" ) ) return( NULL );
 if( FieldPresent( LinePos, "rad" ) && FieldPresent( LinePos, "deg" ) ) return( NULL );
 
 if( FieldPresent( LinePos, "rad" ) )
 {
  if( !GetDouble( LinePos, "rad", angle ) ) return( NULL ); 
 } 
 else
 {
  if( !GetDouble( LinePos, "deg", angle ) ) return( NULL );  
  angle = angle / 180.0 * CONST_PI;
 } 

 if( !FieldPresent( LinePos, "f" ) && !FieldPresent( LinePos, "fcos" ) ) return( NULL );
 if( FieldPresent( LinePos, "f" ) && FieldPresent( LinePos, "fcos" ) ) return( NULL );
 
 if( FieldPresent( LinePos, "f" ) )
 {
  if( !GetDouble( LinePos, "f", factor ) ) return( NULL );
  
  if( fabs( angle ) < CONST_EPSILON )
   factor *= CONST_HUGE_VALUE;
  else
   factor /= ( sin( angle ) * sin( angle ) );
   
 }
 else
  if( !GetDouble( LinePos, "fcos", factor ) ) return( NULL );

 if( factor == 0.0 ) return( NULL );  
   
 return( new HarmonicBendingInteraction( factor, angle, TheData.DefaultParticle.Id, TheData.TheBoxPointer ) );
} 

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

char *StringTheHarmonicBendingInteraction( void *TheInteraction )
{
 HarmonicBendingInteractionPointer TheI = ( HarmonicBendingInteractionPointer )( TheInteraction );
 double f = TheI -> InteractionFactor;
 double angle = acos( TheI -> CosAngle );
 
 sprintf( __ATempStringForSPrintfs, " fcos=%.20le rad=%.20le ", f, angle );
 return( __ATempStringForSPrintfs );
}

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

InteractionPointer ReadThePeriodicTorsionalInteraction( char *LinePos, InputFileDataType& TheData )
{
 double factor, angle;
 int    n = 1;
 
 if( !FieldPresent( LinePos, "rad" ) && !FieldPresent( LinePos, "deg" ) ) return( NULL );
 if( FieldPresent( LinePos, "rad" ) && FieldPresent( LinePos, "deg" ) ) return( NULL );
 
 if( FieldPresent( LinePos, "rad" ) )
 {
  if( !GetDouble( LinePos, "rad", angle ) ) return( NULL ); 
 } 
 else
 {
  if( !GetDouble( LinePos, "deg", angle ) ) return( NULL );  
  angle = angle / 180.0 * CONST_PI;
 } 

 if( !GetDouble( LinePos, "f", factor ) ) return( NULL );
 if( factor == 0.0 ) return( NULL );  
 
 if( FieldPresent( LinePos, "n" ) )
  if( !GetInteger( LinePos, "n", n ) ) return( NULL );
 if( n < 1 ) return( NULL );  
 
 return( new PeriodicTorsionalInteraction( factor, angle, n, TheData.DefaultParticle.Id, TheData.TheBoxPointer ) );
} 

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

char *StringThePeriodicTorsionalInteraction( void *TheInteraction )
{
 PeriodicTorsionalInteractionPointer TheI = ( PeriodicTorsionalInteractionPointer )( TheInteraction );
 double f = TheI -> InteractionFactor;
 double angle = TheI -> Angle;
 int n = TheI -> N;
 
 sprintf( __ATempStringForSPrintfs, " n=%d f=%.20le rad=%.20le ", n, f, angle );
 return( __ATempStringForSPrintfs );
}

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

ConstraintPointer ReadTheDistanceConstraint( char *LinePos, InputFileDataType& TheData )
{
 double distance = -1.0;
 double tolerance = CONST_EPSILON;
 int maxitterations = 1000;
 
 if( FieldPresent( LinePos, "r" ) )
 {
  if( !GetDouble( LinePos, "r", distance ) ) return( NULL );
  if( distance <= 0.0 ) return( NULL );
 } 
 
 if( FieldPresent( LinePos, "error" ) )
 {
  if( !GetDouble( LinePos, "error", tolerance ) ) return( NULL );
  if( tolerance <= 0.0 ) return( NULL );
 } 

 if( FieldPresent( LinePos, "maxstep" ) )
 {
  if( !GetInteger( LinePos, "maxstep", maxitterations ) ) return( NULL );
  if( maxitterations < 1 ) return( NULL );
 } 
  
 return( new DistanceConstraint( distance, tolerance, maxitterations, TheData.DefaultParticle.Id ) );
} 

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

char *StringTheDistanceConstraint( void *TheConstraint )
{
 DistanceConstraintPointer TheC = ( DistanceConstraintPointer )( TheConstraint );
 double distance = TheC -> Distance;
 double tolerance = TheC -> Tolerance;
 int    maxitterations = TheC -> MaxItterations;
 
 sprintf( __ATempStringForSPrintfs, " r=%.20le tol=%.20le maxstep=%u ", distance, tolerance, maxitterations );
 return( __ATempStringForSPrintfs );
}

/* ----------------------------------------------------------------------------- */
/* the table with all known types of interactions and which functions to use    */

typedef InteractionPointer (*__PointerToInteractionReadFunction)( char *LinePos, InputFileDataType& TheData );
typedef char *(*__PointerToInteractionStringFunction)( void *TheInteraction );
typedef struct { InteractionType Type;
                 char *interaction;
                 int length_interaction;
		 __PointerToInteractionReadFunction interaction_read_function;
                 __PointerToInteractionStringFunction interaction_string_function;
	       } __InteractionTableEntry;

__InteractionTableEntry __TheInteractionTable[] = {  { ACoulombInteraction,           "coulomb ",       8,  &ReadTheCoulombInteraction,           &StringTheCoulombInteraction },
                                                     { AGravitationalInteraction,     "gravity ",       8,  &ReadTheGravitationalInteraction,     &StringTheGravitationalInteraction  },
						     { AHarmonicInteraction,          "harmonic ",      9,  &ReadTheHarmonicInteraction,          &StringTheHarmonicInteraction },
						     { ARydbergInteraction,           "rydberg ",       8,  &ReadTheRydbergInteraction,           &StringTheRydbergInteraction  },
						     { AMorseInteraction,             "morse ",         6,  &ReadTheMorseInteraction,             &StringTheMorseInteraction    },
						     { ALennardJonesInteraction,      "lennardjones ", 13,  &ReadTheLennardJonesInteraction,      &StringTheLennardJonesInteraction },
						     { AHarmonicBendingInteraction,   "bending ",       8,  &ReadTheHarmonicBendingInteraction,   &StringTheHarmonicBendingInteraction },
						     { APeriodicTorsionalInteraction, "torsional ",     10, &ReadThePeriodicTorsionalInteraction, &StringThePeriodicTorsionalInteraction }
					          };

#define __TheSizeOfTheInteractionTable (sizeof(__TheInteractionTable)/sizeof(__InteractionTableEntry))

/* ----------------------------------------------------------------------------- */
/* the table with all known types of constraints and which functions to use    */

typedef ConstraintPointer (*__PointerToConstraintReadFunction)( char *LinePos, InputFileDataType& TheData );
typedef char *(*__PointerToConstraintStringFunction)( void *TheConstraint );
typedef struct { ConstraintType Type;
                 char *constraint;
                 int length_constraint;
		 __PointerToConstraintReadFunction constraint_read_function;
                 __PointerToConstraintStringFunction constraint_string_function;
	       } __ConstraintTableEntry;

__ConstraintTableEntry __TheConstraintTable[] = {  { AGeneralDistanceConstraint, "distance ", 9, &ReadTheDistanceConstraint, &StringTheDistanceConstraint }
 				                };

#define __TheSizeOfTheConstraintTable (sizeof(__TheConstraintTable)/sizeof(__ConstraintTableEntry))

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

int ReadTheParticle( char *LinePos, InputFileDataType& TheData )
{
 ParticlePointer NewParticle;
 int ColorIndex, n;

 if( TheData.Periodicity < 0 ) return( 0 );   // if box statement has not been given yet, no particles are allowed...
 if( TheData.DoDynamics ) return ( 0 );
 if( TheData.DoConformation ) return( 0 );
 
 ++TheData.DefaultParticle.Id;                  // by default increase particle ID automatically 

 if( !FieldPresent( LinePos, "x" ) & !FieldPresent( LinePos, "s" ) ) return( 0 );    // one of these fields should always be present
 if( FieldPresent( LinePos, "x" ) & FieldPresent( LinePos, "s" ) ) return( 0 );      // but not both
 if( FieldPresent( LinePos, "p" ) & FieldPresent( LinePos, "ps" ) ) return( 0 );
 
 if( FieldPresent( LinePos, "x" ) )
  if( !GetVector( LinePos, "x", TheData.DefaultParticle.Position ) ) return( 0 );        

 if( FieldPresent( LinePos, "s" ) )
 {
  if( !GetVector( LinePos, "s", TheData.DefaultParticle.Position ) ) return( 0 );         
  TheData.DefaultParticle.Position = SphericalToCartesian( TheData.DefaultParticle.Position.X, TheData.DefaultParticle.Position.Y, TheData.DefaultParticle.Position.Z );
 } 
  
 if( FieldPresent( LinePos, "id" ) )
  if( !GetInteger( LinePos, "id", TheData.DefaultParticle.Id ) ) return( 0 );

 if( FieldPresent( LinePos, "m" ) ) 
  if( !GetDouble( LinePos, "m", TheData.DefaultParticle.Mass ) ) return( 0 );  
 if( TheData.DefaultParticle.Mass <= 0.0 ) return( 0 );
 
 if( FieldPresent( LinePos, "q" ) )  
  if( !GetDouble( LinePos, "q", TheData.DefaultParticle.Charge ) ) return( 0 );                          

 if( FieldPresent( LinePos, "r" ) ) 
  if( !GetDouble( LinePos, "r", TheData.DefaultParticle.Radius ) ) return( 0 );   
 if( TheData.DefaultParticle.Radius <= 0.0 ) return( 0 );                        

 if( FieldPresent( LinePos, "p" ) ) 
  if( !GetVector( LinePos, "p", TheData.DefaultParticle.Momentum ) ) return( 0 );
  
 if( FieldPresent( LinePos, "ps" ) ) 
 {
  if( !GetVector( LinePos, "ps", TheData.DefaultParticle.Momentum ) ) return( 0 );  
  TheData.DefaultParticle.Momentum = SphericalToCartesian( TheData.DefaultParticle.Momentum.X, TheData.DefaultParticle.Momentum.Y, TheData.DefaultParticle.Momentum.Z );
 } 

 if( FieldPresent( LinePos, "c" ) )                                        
  if( GetInteger( LinePos, "c", ColorIndex ) )                  // check if this color is defined
  { 
   if( TheData.Colors == NULL ) return( 0 );
   
   for( n = 0; n < TheData.Colors -> size(); ++n )
    if( (*(TheData.Colors))[ n ] != NULL )
     if( (*(TheData.Colors))[ n ] -> Id == ColorIndex )
      break;
 
   if( n >= TheData.Colors -> size() ) return( 0 ); 
   
   TheData.DefaultParticle.ColorIndex = ColorIndex;   
  }
  else
   return( 0 );
  
 // transform position to rotation, translation and scale
 if( TheData.RotationAngle != 0.0 )
  RotateVectorAroundAxis( TheData.DefaultParticle.Position, TheData.RotationAxis, TheData.RotationAngle );
  
 TheData.DefaultParticle.Position.X += TheData.Translation.X;
 TheData.DefaultParticle.Position.Y += TheData.Translation.Y;
 TheData.DefaultParticle.Position.Z += TheData.Translation.Z;
 
 TheData.DefaultParticle.Position.X *= TheData.Scale.X;
 TheData.DefaultParticle.Position.Y *= TheData.Scale.Y;
 TheData.DefaultParticle.Position.Z *= TheData.Scale.Z;
 
 NewParticle = new Particle( TheData.DefaultParticle );
   
 if( NewParticle == NULL ) return( 0 );

 if( !AddParticleToList( TheData.Particles, (*NewParticle) ) )            // if error or particle already in the list
 {
  delete NewParticle;
  return( 0 );
 }
 
 return( 1 ); 
}

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

int ReadTheInteraction( char *LinePos, InputFileDataType& TheData )
{
 InteractionPointer TheInteractionPointer;
 int n, NrOfParticles, Spot, StopListRead;
 unsigned int ParticleId, FromId, ToId;
 
 if( TheData.Periodicity < 0 ) return( 0 );   // if box statement has not been given yet, no interactions are allowed...
  
 // increment the default id 
 ++TheData.DefaultParticle.Id;
 
 // then check if the id field is present and read it if needed EOS
 if( FieldPresent( LinePos, "id" ) )
  if( !GetInteger( LinePos, "id", TheData.DefaultParticle.Id ) ) return( 0 );
   
 // then try and read the intreaction statement but first remove spaces before the statement
 while( isspace( (int)( *LinePos ) ) ) LinePos++;  
 
 for( n = 0; n < __TheSizeOfTheInteractionTable; ++n )
  if( strstr( LinePos, __TheInteractionTable[ n ].interaction ) == LinePos )                       
  {
   TheInteractionPointer = (*__TheInteractionTable[ n ].interaction_read_function)( LinePos + __TheInteractionTable[ n ].length_interaction, TheData );
   if( TheInteractionPointer == NULL )               
    return( 0 );
   break;                           // break the tables for loop
  }                    
    
 if( n >= __TheSizeOfTheInteractionTable ) return( 0 );                       // check if table entry was found
 
 if( !AddInteractionToList( TheData.Interactions, (*TheInteractionPointer) ) )  // now include it into the list
 {
  delete TheInteractionPointer;                                                 // if it was already included
  return( 0 );
 }
 
 // now start reading the conectivity list, which particles are included in this interaction
 
 LinePos = strstr( LinePos, "{" );                          // locate the { 
 if( LinePos == NULL )
 {
  DeleteInteractionFromList( TheData.Interactions, TheData.DefaultParticle.Id );
  return( 0 );
 } 
 ++LinePos; 
 
 StopListRead = 0;
 Spot = -1;                                                    // set spot to invalid start value
 while( !StopListRead )                                        // continue until we hit } or anything but a number
 {
 
  if( sscanf( LinePos, " %u ", &ParticleId ) == 1 )             // try and read first number
  {                              
  
   // now check if it exists in our data list
   for( Spot = 0; Spot < TheData.Particles -> size(); ++Spot )
    if( (*(TheData.Particles))[ Spot ] != NULL )
     if( (*(TheData.Particles))[ Spot ] -> Id == ParticleId )
     {
      // we have a valid particle ID and the particle is present

      TheInteractionPointer -> IncludeParticle( (*((*(TheData.Particles))[ Spot ])) );
      break;
     }
   
     
   if( Spot >= TheData.Particles -> size() )                 // particle id not defined
   {
    DeleteInteractionFromList( TheData.Interactions, TheData.DefaultParticle.Id );
    return( 0 );
   }
    
   while( isspace( (int)( *LinePos ) ) ) LinePos++;           // now skip all spaces until we hit the first char of the number
   while( isdigit( (int)( *LinePos ) ) ) LinePos++;           // then skip all digits of the number we just read
   while( isspace( (int)( *LinePos ) ) ) LinePos++;           // and again all possible spaces to the next possible number or }
   
   continue;
  }
  else
  {
  
   while( isspace( (int)( *LinePos ) ) ) LinePos++;           // if none found, remove possible spaces
   
   StopListRead = 1;                                          // and we either have hit the } or another non-digit

   if( !strncmp( LinePos, "all", 3 ) )                       // we found the "all" specifier
   {
    for( Spot = 0; Spot < TheData.Particles -> size(); ++Spot )                  // all particles will be included
      if( (*(TheData.Particles))[ Spot ] != NULL )
       TheInteractionPointer -> IncludeParticle( (*((*(TheData.Particles))[ Spot ])) );   

    LinePos += 3;                                               // update LinePos for the next thingy
    while( isspace( (int)( *LinePos ) ) ) LinePos++; 
    StopListRead = 0;
    Spot = -1;                                     // and invalidade Spot
    
    continue;
   }
   
      
   if( !strncmp( LinePos, "to", 2 ) )                        // we found the "to" for a range
   {
    if( Spot <= -1 ) break;                   // we must have read a valid particle ID before
    
    FromId = ParticleId;                                       // the last read ID will be the "from" id
    
    if( sscanf( LinePos + 2, " %u ", &ToId ) != 1 ) break;           // quit if error
    
    if( FromId > ToId ) 
     for( ParticleId = FromId; ParticleId >= ToId; --ParticleId )               // loop over the range of ID's
     {
           
       for( Spot = 0; Spot < TheData.Particles -> size(); ++Spot )                     // and check if it exists in our data list
        if( (*(TheData.Particles))[ Spot ] != NULL )
         if( (*(TheData.Particles))[ Spot ] -> Id == ParticleId )
         {      
          TheInteractionPointer -> IncludeParticle( (*((*(TheData.Particles))[ Spot ])) );    // we have a valid particle ID and the particle is present
          break;
         }
      
       if( Spot >= TheData.Particles -> size() )                 // particle id not defined
       {
        DeleteInteractionFromList( TheData.Interactions, TheData.DefaultParticle.Id );
        return( 0 );
       }	
     }
    else
     for( ParticleId = FromId; ParticleId <= ToId; ++ParticleId )               // loop over the range of ID's
     {
           
       for( Spot = 0; Spot < TheData.Particles -> size(); ++Spot )                     // and check if it exists in our data list
        if( (*(TheData.Particles))[ Spot ] != NULL )
         if( (*(TheData.Particles))[ Spot ] -> Id == ParticleId )
         {      
          TheInteractionPointer -> IncludeParticle( (*((*(TheData.Particles))[ Spot ])) );    // we have a valid particle ID and the particle is present
          break;
         }
      
       if( Spot >= TheData.Particles -> size() )                 // particle id not defined
       {
        DeleteInteractionFromList( TheData.Interactions, TheData.DefaultParticle.Id );
        return( 0 );
       }	
     }
   
    LinePos += 2; 
    while( isspace( (int)( *LinePos ) ) ) LinePos++;                        // else update LinePos to next thingy in the line
    while( isdigit( (int)( *LinePos ) ) ) LinePos++;
    while( isspace( (int)( *LinePos ) ) ) LinePos++;
    StopListRead = 0;
    Spot = -1;                                                              // and invalidade Spot
    
    continue;
   }
   
   
  } 
  
 }
 
 if( (*LinePos) != '}' )                  // we should have stopped at the } char
 {
  DeleteInteractionFromList( TheData.Interactions, TheData.DefaultParticle.Id );
  return( 0 );
 }
 
 return( 1 ); 
}

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

int ReadTheDynamics( char *LinePos, InputFileDataType& TheData )
{
 unsigned int nSteps;
 
 if( TheData.Periodicity < 0 ) return( 0 );   // if box statement has not been given yet, no dynamics yet ...
 if( TheData.DoDynamics ) return( 0 );

 if( !GetDouble( LinePos, "dt", TheData.TimeStep ) ) return( 0 );      // this one should be there
 if( TheData.TimeStep < CONST_EPSILON ) return( 0 );
  
 if( FieldPresent( LinePos, "t" ) )
  if( !GetDouble( LinePos, "t", TheData.Time ) ) return( 0 );
  
 if( FieldPresent( LinePos, "error" ) )
  if( !GetDouble( LinePos, "error", TheData.AccuracyEnergy ) ) return( 0 ); 
 if( TheData.AccuracyEnergy < CONST_EPSILON ) return( 0 );

 if( FieldPresent( LinePos, "n" ) )                                      // either n or tend should be specified
 {
  if( !GetInteger( LinePos, "n", nSteps ) ) return( 0 ); 
  if( nSteps <= 0 ) return( 0 );
  TheData.EndTime = ( double )( nSteps ) * TheData.TimeStep + TheData.Time;
 }
 else
 {
  if( !GetDouble( LinePos, "tend", TheData.EndTime ) ) return ( 0 );
  if( TheData.EndTime < TheData.Time + TheData.TimeStep ) return( 0 );
 }

 TheData.NrOfAutoSnapshotsToTake = 0;
 TheData.AutoSnapShotTimeStep = ( TheData.EndTime - TheData.Time );
 TheData.NextAutoSnapShotTime = TheData.EndTime;
 
 if( FieldPresent( LinePos, "snapshots" ) )
 {
  if( !GetInteger( LinePos, "snapshots", nSteps ) ) return( 0 ); 
  
  if( nSteps < 0 ) return( 0 );
  if( nSteps >= floor( ( TheData.EndTime - TheData.Time ) / TheData.TimeStep ) ) return( 0 );
  
  if( nSteps > 0 )
  {
   TheData.AutoSnapShotTimeStep = ( TheData.EndTime - TheData.Time ) / ( nSteps + 1 );
   TheData.NextAutoSnapShotTime = TheData.Time + TheData.AutoSnapShotTimeStep;
   TheData.NrOfAutoSnapshotsToTake = nSteps;
  }
 }
   
 TheData.DoDynamics = 1; 
 return( 1 ); 
}

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

int ReadTheColor( char *LinePos, InputFileDataType& TheData )
{
 int           Id;
 Vector        RGB;
 ColorPointer  NewColor;
  
 if( !GetInteger( LinePos, "id", Id ) ) return( 0 );                          // these should be here
 if( !GetVector( LinePos, "rgb", RGB ) ) return( 0 );

 NewColor = new Color;
 
 if( NewColor == NULL ) return( 0 );
 
 NewColor -> Id = Id;
 NewColor -> RGB = RGB;
 
 if( !AddColorToList( TheData.Colors, (*NewColor) ) )          // if error, or if color already defined ...
 {
  delete NewColor;
  return( 0 );
 } 
  
 return( 1 ); 
}

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

int ReadTheConformation( char *LinePos, InputFileDataType& TheData )
{
 if( TheData.Periodicity < 0 ) return( 0 );   // if box statement has not been given yet, no steepest descent yet ...
 if( TheData.DoConformation ) return( 0 );
 if( TheData.DoDynamics ) return( 0 );
 
 if( !GetInteger( LinePos, "n", TheData.MaxStepsConformation ) ) return( 0 );               // this should be here
 if( TheData.MaxStepsConformation <= 0 ) return( 0 );
 
 if( FieldPresent( LinePos, "error" ) )
  if( !GetDouble( LinePos, "error", TheData.AccuracyConformation ) ) return( 0 );     
 if( TheData.AccuracyConformation < CONST_EPSILON ) return( 0 );
  
 if( FieldPresent( LinePos, "maxstep" ) )
  if( !GetDouble( LinePos, "maxstep", TheData.MaxStepSizeConformation ) ) return( 0 );  
 if( TheData.MaxStepSizeConformation < CONST_EPSILON ) return( 0 );
 
 TheData.DoConformation = 1;
 return( 1 );
}

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

int ReadTheTranslation( char *LinePos, InputFileDataType& TheData )
{
 if( !FieldPresent( LinePos, "x" ) & !FieldPresent( LinePos, "y" ) & !FieldPresent( LinePos, "z" ) ) return( 0 );
 
 if( FieldPresent( LinePos, "x" ) )
  if( !GetDouble( LinePos, "x", TheData.Translation.X ) ) return( 0 ); 
  
 if( FieldPresent( LinePos, "y" ) )
  if( !GetDouble( LinePos, "y", TheData.Translation.Y ) ) return( 0 );   
  
 if( FieldPresent( LinePos, "z" ) )
  if( !GetDouble( LinePos, "z", TheData.Translation.Z ) ) return( 0 );   
 
 return( 1 ); 
}

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

int ReadTheScale( char *LinePos, InputFileDataType& TheData )
{
 if( !FieldPresent( LinePos, "x" ) & !FieldPresent( LinePos, "y" ) & !FieldPresent( LinePos, "z" ) ) return( 0 );
 
 if( FieldPresent( LinePos, "x" ) )
  if( !GetDouble( LinePos, "x", TheData.Scale.X ) ) return( 0 ); 
  
 if( FieldPresent( LinePos, "y" ) )
  if( !GetDouble( LinePos, "y", TheData.Scale.Y ) ) return( 0 );   
  
 if( FieldPresent( LinePos, "z" ) )
  if( !GetDouble( LinePos, "z", TheData.Scale.Z ) ) return( 0 );   
 
 return( 1 ); 
}

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

int ReadTheRotation( char *LinePos, InputFileDataType& TheData )
{
 if( !FieldPresent( LinePos, "rad" ) & !FieldPresent( LinePos, "deg" ) ) return( 0 );  // one of these should be present
 if( FieldPresent( LinePos, "rad" ) & FieldPresent( LinePos, "deg" ) ) return( 0 );
 
 if( FieldPresent( LinePos, "rad" ) ) 
  if( !GetDouble( LinePos, "rad", TheData.RotationAngle ) ) return( 0 );
 
 if( FieldPresent( LinePos, "deg" ) ) 
 {
  if( !GetDouble( LinePos, "deg", TheData.RotationAngle ) ) return( 0 ); 
  TheData.RotationAngle *= CONST_PI / 180.0;
 }
  
 if( !GetVector( LinePos, "axis", TheData.RotationAxis ) ) return( 0 );
 if( Abs( TheData.RotationAxis ) <= 0.0 ) return( 0 );
 
 return( 1 );
}

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

int ReadTheBox( char *LinePos, InputFileDataType& TheData )
{
 if( TheData.Periodicity >= 0 ) return( 0 );                 // if the box statement has been given already

 if( !FieldPresent( LinePos, "periodic" ) & !FieldPresent( LinePos, "cell" ) ) return( 0 );  // one of these should be present
 if( FieldPresent( LinePos, "periodic" ) & FieldPresent( LinePos, "cell" ) ) return( 0 );
 
 if( FieldPresent( LinePos, "periodic" ) ) 
 {
  if( !GetVector( LinePos, "periodic", TheData.TheBox ) ) return( 0 );
  TheData.Periodicity = 1;
  TheData.TheBoxPointer = &TheData.TheBox;
 } 
 else
 {
  if( !GetVector( LinePos, "cell", TheData.TheBox ) ) return( 0 );
  TheData.Periodicity = 0;
  TheData.TheBoxPointer = NULL;
 }
 
 if( TheData.TheBox.X < CONST_EPSILON ) return( 0 );
 if( TheData.TheBox.Y < CONST_EPSILON ) return( 0 );
 if( TheData.TheBox.Z < CONST_EPSILON ) return( 0 );
 
 return( 1 );
}

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

int ReadTheTemperature( char *LinePos, InputFileDataType& TheData )
{
 if( TheData.DoDynamics ) return ( 0 );
 
 if( !FieldPresent( LinePos, "k" ) ) return( 0 );
 if( !FieldPresent( LinePos, "initial" ) & !FieldPresent( LinePos, "constant" ) ) return( 0 );
 if( FieldPresent( LinePos, "initial" ) & FieldPresent( LinePos, "tau" ) ) return( 0 );
 if( FieldPresent( LinePos, "initial" ) & FieldPresent( LinePos, "rmf" ) ) return( 0 );
 if( FieldPresent( LinePos, "initial" ) & FieldPresent( LinePos, "gamma" ) ) return( 0 );
 if( FieldPresent( LinePos, "initial" ) & FieldPresent( LinePos, "cmrf" ) ) return( 0 );
 if( FieldPresent( LinePos, "rmf" ) & FieldPresent( LinePos, "gamma" ) ) return( 0 );
  
 if( !GetDouble( LinePos, "k", TheData.BoltzmannConstant ) ) return( 0 ); 
 if( TheData.BoltzmannConstant < CONST_EPSILON ) return( 0 );
 
 if( FieldPresent( LinePos, "initial" ) )
 {
  if( !GetDouble( LinePos, "initial", TheData.TargetTemperature ) ) return( 0 );
  if( TheData.TargetTemperature < CONST_EPSILON ) return( 0 );
  
  TheData.ConstantTemperature = -1;
 }
 else
 {
  if( !GetDouble( LinePos, "constant", TheData.TargetTemperature ) ) return( 0 );
  if( TheData.TargetTemperature < CONST_EPSILON ) return( 0 );
  
  TheData.ConstantTemperature |= 1;

  if( FieldPresent( LinePos, "gamma" ) )
  {
   if( !GetDouble( LinePos, "gamma", TheData.GammaFactor ) ) return( 0 );
   if( TheData.GammaFactor < CONST_EPSILON ) return( 0 );
   // if( TheData.GammaFactor > 1.0 ) return( 0 );
   TheData.ConstantTemperature |= 16;
  }

  if( FieldPresent( LinePos, "cmrf" ) ) 
  {
   if( !GetDouble( LinePos, "cmrf", TheData.CenterOfMassRemovalFactor ) ) return( 0 );
   if( TheData.CenterOfMassRemovalFactor < CONST_EPSILON ) return( 0 );
   if( TheData.CenterOfMassRemovalFactor > 1.0 ) return( 0 );
   TheData.ConstantTemperature |= 8;
  }

  if( FieldPresent( LinePos, "rmf" ) ) 
  {
   if( !GetDouble( LinePos, "rmf", TheData.RandomMomentaFactor ) ) return( 0 );
   if( TheData.RandomMomentaFactor < CONST_EPSILON ) return( 0 );
   TheData.ConstantTemperature |= 4;
  }

  if( FieldPresent( LinePos, "tau" ) )
  {
   if( !GetDouble( LinePos, "tau", TheData.TemperatureReScaleTime ) ) return( 0 );
   if( TheData.TemperatureReScaleTime < CONST_EPSILON ) return( 0 );
   TheData.ConstantTemperature |= 2;
  }

 }
 
 return( 1 );
}

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

int ReadTheConstraint( char *LinePos, InputFileDataType& TheData )
{
 ConstraintPointer TheConstraintPointer;
 int n, NrOfParticles, Spot, StopListRead;
 unsigned int ParticleId, FromId, ToId;
 
 if( TheData.Periodicity < 0 ) return( 0 );   // if box statement has not been given yet, no constraints are allowed...
  
 // increment the default id 
 ++TheData.DefaultParticle.Id;
 
 // then check if the id field is present and read it if needed EOS
 if( FieldPresent( LinePos, "id" ) )
  if( !GetInteger( LinePos, "id", TheData.DefaultParticle.Id ) ) return( 0 );
   
 // then try and read the constraint statement but first remove spaces before the statement
 while( isspace( (int)( *LinePos ) ) ) LinePos++;  
 
 for( n = 0; n < __TheSizeOfTheConstraintTable; ++n )
  if( strstr( LinePos, __TheConstraintTable[ n ].constraint ) == LinePos )                       
  {
   TheConstraintPointer = (*__TheConstraintTable[ n ].constraint_read_function)( LinePos + __TheConstraintTable[ n ].length_constraint, TheData );
   if( TheConstraintPointer == NULL )               
    return( 0 );
   break;                           // break the tables for loop
  }                    
    
 if( n >= __TheSizeOfTheConstraintTable ) return( 0 );                       // check if table entry was found
 
 if( !AddConstraintToList( TheData.Constraints, (*TheConstraintPointer) ) )  // now include it into the list
 {
  delete TheConstraintPointer;                                                 // if it was already included
  return( 0 );
 }
 
 // now start reading the conectivity list, which particles are included in this constraint
 
 LinePos = strstr( LinePos, "{" );                          // locate the { 
 if( LinePos == NULL )
 {
  DeleteConstraintFromList( TheData.Constraints, TheData.DefaultParticle.Id );
  return( 0 );
 } 
 ++LinePos; 
 
 StopListRead = 0;
 Spot = -1;                                                    // set spot to invalid start value
 while( !StopListRead )                                        // continue until we hit } or anything but a number
 {
 
  if( sscanf( LinePos, " %u ", &ParticleId ) == 1 )             // try and read first number
  {                              
  
   // now check if it exists in our data list
   for( Spot = 0; Spot < TheData.Particles -> size(); ++Spot )
    if( (*(TheData.Particles))[ Spot ] != NULL )
     if( (*(TheData.Particles))[ Spot ] -> Id == ParticleId )
     {
      // we have a valid particle ID and the particle is present

      TheConstraintPointer -> IncludeParticle( (*((*(TheData.Particles))[ Spot ])) );
      break;
     }
   
     
   if( Spot >= TheData.Particles -> size() )                 // particle id not defined
   {
    DeleteConstraintFromList( TheData.Constraints, TheData.DefaultParticle.Id );
    return( 0 );
   }
    
   while( isspace( (int)( *LinePos ) ) ) LinePos++;           // now skip all spaces until we hit the first char of the number
   while( isdigit( (int)( *LinePos ) ) ) LinePos++;           // then skip all digits of the number we just read
   while( isspace( (int)( *LinePos ) ) ) LinePos++;           // and again all possible spaces to the next possible number or }
   
   continue;
  }
  else
  {
  
   while( isspace( (int)( *LinePos ) ) ) LinePos++;           // if none found, remove possible spaces
   
   StopListRead = 1;                                          // and we either have hit the } or another non-digit
 
   if( !strncmp( LinePos, "to", 2 ) )                        // we found the "to" for a range
   {
    if( Spot <= -1 ) break;                   // we must have read a valid particle ID before
    
    FromId = ParticleId;                                       // the last read ID will be the "from" id
    
    if( sscanf( LinePos + 2, " %u ", &ToId ) != 1 ) break;           // quit if error
    
    if( FromId > ToId ) 
     for( ParticleId = FromId; ParticleId >= ToId; --ParticleId )               // loop over the range of ID's
     {
           
       for( Spot = 0; Spot < TheData.Particles -> size(); ++Spot )                     // and check if it exists in our data list
        if( (*(TheData.Particles))[ Spot ] != NULL )
         if( (*(TheData.Particles))[ Spot ] -> Id == ParticleId )
         {      
          TheConstraintPointer -> IncludeParticle( (*((*(TheData.Particles))[ Spot ])) );    // we have a valid particle ID and the particle is present
          break;
         }
      
       if( Spot >= TheData.Particles -> size() )                 // particle id not defined
       {
        DeleteConstraintFromList( TheData.Constraints, TheData.DefaultParticle.Id );
        return( 0 );
       }	
     }
    else
     for( ParticleId = FromId; ParticleId <= ToId; ++ParticleId )               // loop over the range of ID's
     {
           
       for( Spot = 0; Spot < TheData.Particles -> size(); ++Spot )                     // and check if it exists in our data list
        if( (*(TheData.Particles))[ Spot ] != NULL )
         if( (*(TheData.Particles))[ Spot ] -> Id == ParticleId )
         {      
          TheConstraintPointer -> IncludeParticle( (*((*(TheData.Particles))[ Spot ])) );    // we have a valid particle ID and the particle is present
          break;
         }
      
       if( Spot >= TheData.Particles -> size() )                 // particle id not defined
       {
        DeleteConstraintFromList( TheData.Constraints, TheData.DefaultParticle.Id );
        return( 0 );
       }	
     }
   
    LinePos += 2; 
    while( isspace( (int)( *LinePos ) ) ) LinePos++;                        // else update LinePos to next thingy in the line
    while( isdigit( (int)( *LinePos ) ) ) LinePos++;
    while( isspace( (int)( *LinePos ) ) ) LinePos++;
    StopListRead = 0;
    Spot = -1;                                                              // and invalidade Spot
    
    continue;
   }
   
   
  } 
  
 }
 
 if( (*LinePos) != '}' )                  // we should have stopped at the } char
 {
  DeleteConstraintFromList( TheData.Constraints, TheData.DefaultParticle.Id );
  return( 0 );
 }
 
 return( 1 ); 
}

/* ----------------------------------------------------------------------------- */
/*   the table with all known types of statements and which functions to use     */

typedef int (*__PointerToStatementFunction)( char *LinePos, InputFileDataType& TheData );
typedef struct { char *statement; int length_statement; __PointerToStatementFunction statement_read_function; } __StatementTableEntry;

__StatementTableEntry __TheStatementTable[] = {  { "particle ",      9,  &ReadTheParticle },
                                                 { "interaction ",  12,  &ReadTheInteraction },
					         { "color ",         6,  &ReadTheColor },
					         { "dynamics ",      9,  &ReadTheDynamics },
						 { "conformation ", 13,  &ReadTheConformation },
						 { "translate ",    10,  &ReadTheTranslation },
						 { "scale ",         5,  &ReadTheScale },      
						 { "rotate ",        7,  &ReadTheRotation },
						 { "box ",           4,  &ReadTheBox },
						 { "temperature ",  12,  &ReadTheTemperature },
						 { "constrain ",    10,  &ReadTheConstraint },
					      };

#define __TheSizeOfTheStatementTable (sizeof(__TheStatementTable)/sizeof(__StatementTableEntry))

/* ---------------------------------------- The public available routines ------ */

int DeleteInteractionFromList( InteractionListPointer& Interactions, unsigned int Id )
{
 InteractionPointer TempPointer;
 
 if( Interactions == NULL ) return( 0 );
 
 for( int n = 0; n < Interactions -> size(); ++n )
  if( (*(Interactions))[ n ] != NULL )
   if( (*(Interactions))[ n ] -> Id == Id )
   {
    TempPointer = (*(Interactions))[ n ];
    (*(Interactions))[ n ] = NULL;
    delete TempPointer;
    return( 1 );
   }
  
 return( 0 );
}

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

int DeleteParticleFromList( ParticleListPointer& Particles, unsigned int Id )
{
 ParticlePointer TempPointer;
 
 if( Particles == NULL ) return( 0 );
 
 for( int n = 0; n < Particles -> size(); ++n )
  if( (*(Particles))[ n ] != NULL )
   if( (*(Particles))[ n ] -> Id == Id )
   {
    TempPointer = (*(Particles))[ n ];
    (*(Particles))[ n ] = NULL;
    delete TempPointer;
    return( 1 );
   }
  
 return( 0 );
}

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

int DeleteConstraintFromList( ConstraintListPointer& Constraints, unsigned int Id )
{
 ConstraintPointer TempPointer;
 
 if( Constraints == NULL ) return( 0 );
 
 for( int n = 0; n < Constraints -> size(); ++n )
  if( (*(Constraints))[ n ] != NULL )
   if( (*(Constraints))[ n ] -> Id == Id )
   {
    TempPointer = (*(Constraints))[ n ];
    (*(Constraints))[ n ] = NULL;
    delete TempPointer;
    return( 1 );
   }
  
 return( 0 );
}

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

int DeleteColorFromList( ColorListPointer& Colors, int Id )
{
 ColorPointer TempPointer;
 
 if( Colors == NULL ) return( 0 );
 
 for( int n = 0; n < Colors -> size(); ++n )
  if( (*(Colors))[ n ] != NULL )
   if( (*(Colors))[ n ] -> Id == Id )
   {
    TempPointer = (*(Colors))[ n ];
    (*(Colors))[ n ] = NULL;
    delete TempPointer;
    return( 1 );
   }

 return( 0 );
}

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

void DeleteColorList( ColorListPointer& Colors )
{
 if( Colors != NULL ) 
 {
  for( int n = 0; n < Colors -> size(); ++n )
   if( (*(Colors))[ n ] != NULL )
    delete (*(Colors))[ n ];
  delete Colors;
  Colors = NULL;
 }
}

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

void DeleteTheData( InputFileDataType& TheData )
{
 DeleteInteractionList( TheData.Interactions );
 DeleteParticleList( TheData.Particles );
 DeleteColorList( TheData.Colors );
 DeleteConstraintList( TheData.Constraints );
}

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

void DeleteInteractionList( InteractionListPointer& Interactions )
{
 if( Interactions != NULL ) 
 {
  for( int n = 0; n < Interactions -> size(); ++n )
   if( (*(Interactions))[ n ] != NULL )
    delete (*(Interactions))[ n ];
  delete Interactions;
  Interactions = NULL;
 }
}

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

void DeleteParticleList( ParticleListPointer& Particles )
{
 if( Particles != NULL ) 
 {
  for( int n = 0; n < Particles -> size(); ++n )
   if( (*(Particles))[ n ] != NULL )
    delete (*(Particles))[ n ];
  delete Particles;  
  Particles = NULL;
 }
}

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

void DeleteConstraintList( ConstraintListPointer& Constraints )
{
 if( Constraints != NULL ) 
 {
  for( int n = 0; n < Constraints -> size(); ++n )
   if( (*(Constraints))[ n ] != NULL )
    delete (*(Constraints))[ n ];
  delete Constraints;  
  Constraints = NULL;
 }
}

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

int AddColorToList( ColorListPointer& Colors, Color& C )
{
 int n, EmptySpot;
 ColorPointer NullPointer = NULL;
 
 if( Colors == NULL )
 {
  Colors = new ColorList( 5 );
  if( Colors == NULL ) return( 0 ); 
  for( n = 0; n < Colors -> size(); ++n )
   (*(Colors))[ n ] = NULL;
 } 
  
 for( EmptySpot = -1, n = 0; n < Colors -> size(); n++ )
  if( (*(Colors))[ n ] == NULL )
  {
   if( EmptySpot <= -1 )
    EmptySpot = n;
  } 
  else
   if( (*(Colors))[ n ] -> Id == C.Id )
    return( 0 );
     
 if( EmptySpot <= -1 )
 {
  Colors -> insert( Colors -> end(), 5, NullPointer );
  EmptySpot = n;
 }
 
 (*(Colors))[ EmptySpot ] = &C;
 return( 1 );   
}

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

int AddParticleToList( ParticleListPointer& Particles, Particle& P )
{
 int n, EmptySpot;
 ParticlePointer NullPointer = NULL;
 
 if( Particles == NULL )
 {
  Particles = new ParticleList( 20 );
  if( Particles == NULL ) return( 0 ); 
  for( n = 0; n < Particles -> size(); ++n )
   (*(Particles))[ n ] = NULL;
 } 
  
 for( EmptySpot = -1, n = 0; n < Particles -> size(); n++ )
  if( (*(Particles))[ n ] == NULL )
  {
   if( EmptySpot <= -1 )
    EmptySpot = n;
  } 
  else
   if( (*(Particles))[ n ] -> Id == P.Id )
    return( 0 );
     
 if( EmptySpot <= -1 )
 {
  Particles -> insert( Particles -> end(), 5, NullPointer );
  EmptySpot = n;
 }
 
 (*(Particles))[ EmptySpot ] = &P;
 return( 1 );   
}

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

int AddInteractionToList( InteractionListPointer& Interactions, Interaction& I )
{
 int n, EmptySpot;
 InteractionPointer NullPointer = NULL;
 
 if( Interactions == NULL )
 {
  Interactions = new InteractionList( 20 );
  if( Interactions == NULL ) return( 0 ); 
  for( n = 0; n < Interactions -> size(); ++n )
   (*(Interactions))[ n ] = NULL;
 } 
 
 for( EmptySpot = -1, n = 0; n < Interactions -> size(); n++ )
  if( (*(Interactions))[ n ] == NULL )
  {
   if( EmptySpot <= -1 )
    EmptySpot = n;
  }  
  else
   if( (*(Interactions))[ n ] -> Id == I.Id )
    return( 0 );
    
 if( EmptySpot <= -1 )
 {
  Interactions -> insert( Interactions -> end(), 5, NullPointer );
  EmptySpot = n;
 }
 
 (*(Interactions))[ EmptySpot ] = &I;
 return( 1 );   
}

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

int AddConstraintToList( ConstraintListPointer& Constraints, Constraint& C )
{
 int n, EmptySpot;
 ConstraintPointer NullPointer = NULL;
 
 if( Constraints == NULL )
 {
  Constraints = new ConstraintList( 20 );
  if( Constraints == NULL ) return( 0 ); 
  for( n = 0; n < Constraints -> size(); ++n )
   (*(Constraints))[ n ] = NULL;
 } 
 
 for( EmptySpot = -1, n = 0; n < Constraints -> size(); n++ )
  if( (*(Constraints))[ n ] == NULL )
  {
   if( EmptySpot <= -1 )
    EmptySpot = n;
  }  
  else
   if( (*(Constraints))[ n ] -> Id == C.Id )
    return( 0 );
    
 if( EmptySpot <= -1 )
 {
  Constraints -> insert( Constraints -> end(), 5, NullPointer );
  EmptySpot = n;
 }
 
 (*(Constraints))[ EmptySpot ] = &C;
 return( 1 );   
}

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

void SortInteractionList( InteractionListPointer& Interactions )
{
 int n, m, N, Index;
 InteractionPointer TempPointer;

 if( Interactions == NULL ) return;

 for( N = n = 0; n < Interactions -> size(); n++ )     // first count number of items present
  if( (*(Interactions))[ n ] != NULL ) N++;

 for( n = 0; n < N; n++ )
 {
  for( Index = n; Index < Interactions -> size(); Index++ )    // get first non-empty slot 
   if( (*(Interactions))[ Index ] != NULL ) break;

  for( m = Index + 1; m < Interactions -> size(); m++ )       // loop through the rest and find lowest id
   if( (*(Interactions))[ m ] != NULL )
    if( (*(Interactions))[ m ] -> Id < (*(Interactions))[ Index ] -> Id )
     Index = m;

  TempPointer = (*(Interactions))[ n ];                       // swap lowest found value into slot
  (*(Interactions))[ n ] = (*(Interactions))[ Index ];
  (*(Interactions))[ Index ] = TempPointer;
 }

}

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

void SortParticleList( ParticleListPointer& Particles )
{
 int n, m, N, Index;
 ParticlePointer TempPointer;

 if( Particles == NULL ) return;

 for( N = n = 0; n < Particles -> size(); n++ )     // first count number of items present
  if( (*(Particles))[ n ] != NULL ) N++;

 for( n = 0; n < N; n++ )
 {
  for( Index = n; Index < Particles -> size(); Index++ )    // get first non-empty slot
   if( (*(Particles))[ Index ] != NULL ) break;

  for( m = Index + 1; m < Particles -> size(); m++ )       // loop through the rest and find lowest id
   if( (*(Particles))[ m ] != NULL )
    if( (*(Particles))[ m ] -> Id < (*(Particles))[ Index ] -> Id )
     Index = m;

  TempPointer = (*(Particles))[ n ];                       // swap lowest found value into slot
  (*(Particles))[ n ] = (*(Particles))[ Index ];
  (*(Particles))[ Index ] = TempPointer;
 }

}

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

void SortColorList( ColorListPointer& Colors )
{
 int n, m, N, Index;
 ColorPointer TempPointer;

 if( Colors == NULL ) return;

 for( N = n = 0; n < Colors -> size(); n++ )     // first count number of items present
  if( (*(Colors))[ n ] != NULL ) N++;

 for( n = 0; n < N; n++ )
 {
  for( Index = n; Index < Colors -> size(); Index++ )    // get first non-empty slot
   if( (*(Colors))[ Index ] != NULL ) break;

  for( m = Index + 1; m < Colors -> size(); m++ )       // loop through the rest and find lowest id
   if( (*(Colors))[ m ] != NULL )
    if( (*(Colors))[ m ] -> Id < (*(Colors))[ Index ] -> Id )
     Index = m;

  TempPointer = (*(Colors))[ n ];                       // swap lowest found value into slot
  (*(Colors))[ n ] = (*(Colors))[ Index ];
  (*(Colors))[ Index ] = TempPointer;
 }

}

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

void SortConstraintList( ConstraintListPointer& Constraints )
{
 int n, m, N, Index;
 ConstraintPointer TempPointer;

 if( Constraints == NULL ) return;

 for( N = n = 0; n < Constraints -> size(); n++ )     // first count number of items present
  if( (*(Constraints))[ n ] != NULL ) N++;

 for( n = 0; n < N; n++ )
 {
  for( Index = n; Index < Constraints -> size(); Index++ )    // get first non-empty slot
   if( (*(Constraints))[ Index ] != NULL ) break;

  for( m = Index + 1; m < Constraints -> size(); m++ )       // loop through the rest and find lowest id
   if( (*(Constraints))[ m ] != NULL )
    if( (*(Constraints))[ m ] -> Id < (*(Constraints))[ Index ] -> Id )
     Index = m;

  TempPointer = (*(Constraints))[ n ];                       // swap lowest found value into slot
  (*(Constraints))[ n ] = (*(Constraints))[ Index ];
  (*(Constraints))[ Index ] = TempPointer;
 }

}

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

int ReadDataFromFile( char *FileName, InputFileDataType& TheData, int &LineNumber, FILE *InFile )
{
 char TheLineRead[ MAX_LINE_LENGTH ];
 FILE *InputFile = NULL;
 char *LinePos = NULL;
 InteractionPointer NewInteraction;  
 Particle TempParticle;
 Interaction TempInteraction;
 int n;
 
 if( FileName == NULL ) return( 0 );
  
 if( InFile == NULL )
  InputFile = fopen( FileName, "rt" );
 else
  InputFile = InFile; 
  
 if( InputFile == NULL ) return( 0 );
 
 TheData.Periodicity = -1;                // give it an invalid state to detect if box statment has been given before
 TheData.TheBoxPointer = NULL;            // particle or interaction statements are encountered...
  
 LineNumber = 0;
 while( !feof( InputFile ) )
 {
  ++LineNumber;
  
  LinePos = fgets( TheLineRead, sizeof( TheLineRead ), InputFile );            // get a line

  if( LinePos == NULL )
  {
   if( feof( InputFile ) )
    break;
   if( InFile == NULL )
    fclose( InputFile );
   return( 0 ); 
  }
  
  for( n = 0; TheLineRead[ n ] != '\0'; ++n )                       // make it all lowercase
   TheLineRead[ n ] = ( char )( tolower( TheLineRead[ n ] ) );
   
  while( isspace( (int)( *LinePos ) ) ) LinePos++;                              // ignore spaces at beginning
  
  if( LinePos[ 0 ] == '\0' ) continue;                                          // and also ignore empty lines
  if( LinePos[ 0 ] == '\n' ) continue;                        
  if( LinePos[ 0 ] == '\r' ) continue;     

  if( strstr( LinePos, "# end of data" ) != NULL )                    // if we hit the end of data label
  {
   if( InFile == NULL )
    fclose( InputFile );
   return( 1 );
  }
  
  if( LinePos[ 0 ] == '#' ) continue;                                           // ignore comments      
  
  if( strstr( LinePos, "#" ) != NULL ) *strstr( LinePos, "#" ) = '\0';         // remove comments within a line
  
  for( n = 0; n < __TheSizeOfTheStatementTable; ++n )
   if( strstr( LinePos, __TheStatementTable[ n ].statement ) == LinePos )                       
   {
    
    if( !(*__TheStatementTable[ n ].statement_read_function)( LinePos + __TheStatementTable[ n ].length_statement, TheData ) )               
    {
     if( InFile == NULL )
      fclose( InputFile );
     return( 0 );
    }
    
    break;                           // break the tables for loop
   }                    
    
  if( n < __TheSizeOfTheStatementTable ) continue;           // check if table entry was found
    
  // we found an illegal start statement or some error has occured
  if( InFile == NULL )
   fclose( InputFile );
  return( 0 );
 }
 
 if( InFile == NULL )
  fclose( InputFile );
 return( 1 );
}

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

int WriteDataToFile( char *FileName, InputFileDataType& TheData, FILE *OutFile )
{
 char TheOptions[ 512 ];
 ParticlePointer    TheParticle;
 InteractionPointer TheInteraction;
 ConstraintPointer  TheConstraint;
 FILE *OutputFile = NULL;
 int n, Index;

 if( FileName == NULL ) return( 0 );
 
 if( OutFile == NULL )
  OutputFile = fopen( FileName, "a" );
 else
  OutputFile = OutFile;
  
 if( OutputFile == NULL ) return ( 0 );
 
 
 // first output color table if exist
 if( TheData.Colors != NULL )
 {
  SortColorList( TheData.Colors );                  // sort list first

  fputs( "# the color table\n", OutputFile );
  
  for( n = 0; n < TheData.Colors -> size(); ++n )
   if( (*(TheData.Colors))[ n ] != NULL )
    fprintf( OutputFile, "color id=%u rgb=[%.20le,%.20le,%.20le]\n", (*(TheData.Colors))[ n ] -> Id, 
             (*(TheData.Colors))[ n ] -> RGB.X, (*(TheData.Colors))[ n ] -> RGB.Y, (*(TheData.Colors))[ n ] -> RGB.Z ); 
    
  fputs( "\n", OutputFile );
 }
 
 
 // dump the box statement and specify if we used periodicity
 fputs( "# the box statement\n", OutputFile );
 if( TheData.Periodicity )
  fprintf( OutputFile, "box periodic=[%.20le,%.20le,%.20le]\n\n", TheData.TheBox.X, TheData.TheBox.Y, TheData.TheBox.Z );
 else 
  fprintf( OutputFile, "box cell=[%.20le,%.20le,%.20le]\n\n", TheData.TheBox.X, TheData.TheBox.Y, TheData.TheBox.Z );
  
 // then dump all the particles
 if( TheData.Particles != NULL )
 {
  SortParticleList( TheData.Particles );       // sort list first

  fputs( "# all the particles\n", OutputFile );
  
  for( n = 0; n < TheData.Particles -> size(); ++n )
   if( (*(TheData.Particles))[ n ] != NULL )
   {
    TheParticle = (*(TheData.Particles))[ n ];
    fprintf( OutputFile, "particle r=%.20le m=%.20le q=%.20le", TheParticle -> Radius, TheParticle -> Mass, TheParticle -> Charge );
    fprintf( OutputFile, " x=[%.20le,%.20le,%.20le]", TheParticle -> Position.X, TheParticle -> Position.Y, TheParticle -> Position.Z ); 
    fprintf( OutputFile, " p=[%.20le,%.20le,%.20le]", TheParticle -> Momentum.X, TheParticle -> Momentum.Y, TheParticle -> Momentum.Z ); 
    if( TheData.Colors != NULL )
     fprintf( OutputFile, " c=%d", TheParticle -> ColorIndex );
    fprintf( OutputFile, " id=%u\n", TheParticle -> Id );
   }
  
  fputs( "\n", OutputFile );       
 }
 
 
 // and dump the interactions 
 if( TheData.Interactions != NULL )
 {
  SortInteractionList( TheData.Interactions );          // sort list first

  fputs( "# all the interactions\n", OutputFile );
   
  for( n = 0; n < TheData.Interactions -> size(); ++n )
   if( (*(TheData.Interactions))[ n ] != NULL )
   {
    TheInteraction = (*(TheData.Interactions))[ n ];
    
    // lookup in table

    for( Index = 0; Index < __TheSizeOfTheInteractionTable; ++Index )
     if( TheInteraction -> Type == __TheInteractionTable[ Index ].Type )
      break;
    
    // if not known, then stop  
    if( Index >= __TheSizeOfTheInteractionTable )
    {
     if( OutFile == NULL )
      fclose( OutputFile );
     return( 0 );
    }
    
    // and start writing the statement to file
    fprintf( OutputFile, "interaction %sid=%u  %s { ", __TheInteractionTable[ Index ].interaction, TheInteraction -> Id,
                      (*__TheInteractionTable[ Index ].interaction_string_function)( TheInteraction ) ); 
 
    // now dump all it's particles 
    for( Index = 0; Index < TheInteraction -> GetParticleListSize(); ++Index )
     if( TheInteraction -> GetParticlePointer( Index ) != NULL )
      fprintf( OutputFile, "%u ", TheInteraction -> GetParticlePointer( Index ) -> Id );
          
    fprintf( OutputFile, "}\n" );   
   }
  
  fputs( "\n", OutputFile );       
 }
 
 
 // also dump the constraints
 if( TheData.Constraints != NULL )
 {
  SortConstraintList( TheData.Constraints );           // sort list first

  fputs( "# all the constraints\n", OutputFile );
   
  for( n = 0; n < TheData.Constraints -> size(); ++n )
   if( (*(TheData.Constraints))[ n ] != NULL )
   {
    TheConstraint = (*(TheData.Constraints))[ n ];
    
    // lookup in table
    for( Index = 0; Index < __TheSizeOfTheConstraintTable; ++Index )
     if( TheConstraint -> Type == __TheConstraintTable[ Index ].Type )
      break;
    
    // if not known, then stop  
    if( Index >= __TheSizeOfTheConstraintTable )
    {
     if( OutFile == NULL )
      fclose( OutputFile );
     return( 0 );
    }
    
    // and start writing the statement to file
    fprintf( OutputFile, "constrain %sid=%u  %s { ", __TheConstraintTable[ Index ].constraint, TheConstraint -> Id,
                      (*__TheConstraintTable[ Index ].constraint_string_function)( TheConstraint ) ); 
 
    // now dump all it's particles 
    for( Index = 0; Index < TheConstraint -> GetParticleListSize(); ++Index )
     if( TheConstraint -> GetParticleId( Index ) != 0 )
      fprintf( OutputFile, "%u ", TheConstraint -> GetParticleId( Index ) );
          
    fprintf( OutputFile, "}\n" );   
   }
  
  fputs( "\n", OutputFile );       
 }
 
 
 // write conformation stuff
 if( TheData.DoConformation )
 {
  if( TheData.DoConformation < 0 )
   fprintf( OutputFile, "# conformation search statement\n# conformation n=%u maxstep=%.20le error=%.20le\n\n", 
            TheData.MaxStepsConformation, TheData.MaxStepSizeConformation, TheData.AccuracyConformation );
	    
  if( TheData.DoConformation > 0 )
   fprintf( OutputFile, "# conformation search statement\nconformation n=%u maxstep=%.20le error=%.20le\n\n", 
            TheData.MaxStepsConformation, TheData.MaxStepSizeConformation, TheData.AccuracyConformation );
 }
   
 
 
 // write the temperature statement, if needed ...
 if( TheData.ConstantTemperature )
 {
  if( TheData.ConstantTemperature < 0 )    // only init at t=0.0 but use energy conservation
   fprintf( OutputFile, "# if needed, setup the initial temperature\ntemperature k=%.20le initial=%.20le\n\n", 
            TheData.BoltzmannConstant, TheData.TargetTemperature );

  if( TheData.ConstantTemperature > 0 )
  {
   TheOptions[ 0 ] = 0;

   if( TheData.RandomMomentaFactor >= CONST_EPSILON )
    sprintf( TheOptions, "%s rmf=%.20le", TheOptions, TheData.RandomMomentaFactor ); 

   if( TheData.CenterOfMassRemovalFactor >= CONST_EPSILON )
    sprintf( TheOptions, "%s cmrf=%.20le", TheOptions, TheData.CenterOfMassRemovalFactor ); 

   if( TheData.GammaFactor >= CONST_EPSILON )
    sprintf( TheOptions, "%s gamma=%.20le", TheOptions, TheData.GammaFactor ); 

   if( TheData.TemperatureReScaleTime >= CONST_EPSILON )
    sprintf( TheOptions, "%s tau=%.20le", TheOptions, TheData.TemperatureReScaleTime ); 

   fprintf( OutputFile, "# keep the temperature constant\ntemperature k=%.20le constant=%.20le%s\n\n", 
             TheData.BoltzmannConstant, TheData.TargetTemperature, TheOptions ); 
  
  }
 }
 
 // finally write dynamics 
 if( TheData.DoDynamics )
 {
  if( TheData.DoDynamics < 0 )
   fprintf( OutputFile, "# dynamics statement\n# dynamics t=%.20le dt=%.20le tend=%.20le error=%.20le snapshots=%u\n\n",
            TheData.Time, TheData.TimeStep, TheData.EndTime, TheData.AccuracyEnergy, TheData.NrOfAutoSnapshotsToTake );
	    
  if( TheData.DoDynamics > 0 )
   fprintf( OutputFile, "# dynamics statement\ndynamics t=%.20le dt=%.20le tend=%.20le error=%.20le snapshots=%u\n\n",
            TheData.Time, TheData.TimeStep, TheData.EndTime, TheData.AccuracyEnergy, TheData.NrOfAutoSnapshotsToTake );	    
 }

 // and mark end of data
 fprintf( OutputFile, "# end of data\n\n\n" );
 
 if( OutFile == NULL )
  fclose( OutputFile );
 return( 1 );
}

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