/*
   Project: Adun

   Copyright (C) 2005 Michael Johnston & Jordi Villa-Freixa

   Author: Michael Johnston

   Created: 2005-06-23 11:06:55 +0200 by michael johnston

   This application is free software; you can redistribute it and/or
   modify it under the terms of the GNU General Public
   License as published by the Free Software Foundation; either
   version 2 of the License, or (at your option) any later version.

   This application is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Library General Public License for more details.

   You should have received a copy of the GNU General Public
   License along with this library; if not, write to the Free
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111 USA.
*/
#ifndef _ADLANGEVIN_THERMOSTAT_
#define _ADLANGEGIN_THERMOSTAT_

#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "AdunKernel/AdunSimulator.h"

/**
AdLangevinThermostat objects are AdSimulator components that
transform the newtonian dynamics generated by AdSimulator objects into
langevin dynamics.

The velocity verlet algorithm for langevin dynamics is as follows (BBK) - Here
\f$R_{n}\f$ is a random acceleration vector and \f$\gamma\f$ is the collision parameter.

\f[
\vec v_{i+\frac{1}{2}} =  \vec v_{i} + \frac{t}{2} \left (\vec a_{i} - \gamma \vec v_{i} + \vec R_{i} \right ) \\
\f]
\f[
\vec r_{i+1} = \vec r_{i} + \vec v_{i+ \frac{1}{2}}t \\
\f]
\f[
\vec v_{i+1} =  \vec v_{i + \frac{1}{2}} + \frac{t}{2} \left (\vec a_{i+1} - \gamma \vec v_{i+1} + \vec R_{i+1} \right ) \\
\f]

The third equation defines \f$\vec v_{i +1}\f$ implicitly. However the linear dependency allows solution for 
 \f$\vec v_{i +1}\f$ explicitly.

To achieve the above AdLangevinThermostat instances apply the following modifications to the verlet algorithm
defined by AdSimulator. \f$ v{'}\f$ indicates an intermediate form of the velocity vector. The arrows on the
vectors are omitted for clarity.

\f[
v_{i}^{'} =  v_{i} + \frac{t}{2} \left ( R_{i} - \gamma v_{i} \right ) \\
\f]
\f[
v_{i+\frac{1}{2}} =  v_{i}^{'} + \frac{ a_{i}t}{2} \\
\f]
\f[
 r_{i+1} = r_{i} + v_{i+ \frac{1}{2}}t \\
\f]
\f[
v_{i + 1}^{'} = v_{i + \frac{1}{2}} + \frac{a_{i+1}t}{2} \\
\f]
\f[
v_{i+1} =  \frac{ (2 v_{i + 1}^{'} + R_{i+1} t ) }{2 + \gamma t} \\
\f]

The second, third and fourth steps are those performed by the AdSimulator object. The first step
is computed by the AdLangevinThermostat object on receipt of simulatorWillPerformFirstVelocityUpdateForSystem:()
and the last when simulatorDidPerformSecondVelocityUpdateForSystem:() is sent.

The magnitude of the random acceleration vector applied to each element is generated from a gaussian distribution with mean 0
and variance \f$ \frac{6\gamma k_{b} T}{m \Delta t} \f$. where \f$m\f$ is the mass of the element and \f$ \Delta t \f$ 
is the simulation time step. Note that the associated random force vector has mean 0 and variance
\f$6 \frac{\gamma k_{b}mT}{\Delta t}\f$ as required by the fluctuation-dissipation theorem.
The direction is given by a randomly generated unit vector which has a uniform distribution over the surface of the unit-sphere.

\todo Missing Functionality - Enable ability to set a gamma value for each element
\ingroup Inter
**/

@interface AdLangevinThermostat: NSObject <AdSimulatorComponent>
{
	BOOL removeDOF; //!< Indicates if the translation DOF should be periodically be removed
	int removeDOFInterval; //!< The interval at which the removal should be done.
	double timeStep;
	double gamma;
	double targetTemperature;
	double variance;
	gsl_rng* twister;
	NSMutableDictionary* matrixDict;
	NSMutableDictionary* massesDict;
	AdSystemCollection* systemCollection;
	AdForceFieldCollection* forceFieldCollection;
}
/**
As initWithTargetTemperature:() with a temperature of 300K.
*/
- (id) init;
/**
As initWithTargetTemperature:gamma:() with a gamma of 0.05.
*/
- (id) initWithTargetTemperature: (double) doubleOne;
/**
As initWithTargetTemperature:gamma:seed:() with seed set to 322
*/
- (id) initWithTargetTemperature: (double) doubleOne gamma: (double) doubleTwo;
/**
As initWithTargetTemperature:gamma:seed:resetDOFInterval:() with \e interval set to 1000
*/
- (id) initWithTargetTemperature: (double) doubleOne gamma: (double) doubleTwo seed: (int) anInt;
/**
Designated initialiser.
\param doubleOne The temperature to be maintained by the thermostat (Kelvin).
\param doubleTwo The value of gamma. Should be in \f$fs^{-1}\f$.
\param anInt Seed for random number generation
\param interval Interval at which to remove translation DOF from the systems being simulated.
If this is less than zero no removal is done.
*/
- (id) initWithTargetTemperature: (double) doubleOne 
			   gamma: (double) doubleTwo
			    seed: (int) anInt
		removeDOFInterval: (int) interval;
/**
Value of the friction coefficent used for the langevin dynamics.
*/
- (double) gamma;
/**
Set the value of gamma.
*/
- (void) setGamma: (double) value;
/**
Sets the target temperature for the thermostat to \e aDouble.
If \e aDouble is less than 0 an NSInvalidArgumentException is raised.
*/
- (void) setTargetTemperature: (double) aDouble;
/**
Returns the target temperature
*/
- (double) targetTemperature;
/**
Sets the seed used for initialising the random number generator
which generates the random force matrix. Has the effect of
resetting the random number generator if it had already been
used.*/
- (void) setSeed: (int) anInt;
@end

#endif
