//****************************************************************************

//! \file Tsimmain.hpp
//! \date 24-11-2003

//! \class Tsimmain
//! \brief The main class of the simulation

//*****************************************************************************

//     This file is part of heliumvmc.
//
//     heliumvmc 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.
//
//     heliumvmc 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 General Public License for more details.
//
//     You should have received a copy of the GNU General Public License
//     along with heliumvmc; if not, write to the Free Software
//     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

//*****************************************************************************

#ifndef TSIMMAIN_HPP
#define TSIMMAIN_HPP

#include <vector>
#include <numeric>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <functional>
#include <cmath>
#include <ctime>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

const double channelsize = .05;
const int channels = 100;


using std::vector;

//*****************************************************************************
template <class Arg>
class vAddRand:  public std::unary_function<Arg, void>
//-----------------------------------------------------------------------------
//! A unary function object derived from unary_function to add a random
//! number to a vector
//*****************************************************************************
{
	double bound;//! the boundary for the random number distribution
	gsl_rng *gslrand1;//! the gsl random number generator

public:
	//! \brief A constructor does nothing but to initialize the boundary
	//! \param g1 a gsl random number generator
	//! \param x the boundary
	vAddRand( gsl_rng *g1, const double &x) : gslrand1(g1), bound( x) {};

	//! \brief Another constructor doing nothing (if you want to initialize the boundary later)
	vAddRand() : bound( 0.){};
	
	//! \brief resets the boundary
	//! \param b the new boundary
	void resetBound( double b) { bound = b;}

	~vAddRand() {};

	//! \brief resets the boundary
	//! \param x the new boundary
	void mreset( const double &x) { bound = x;}

	//! \brief the function doing the real work
	//! \param x the value to be added to the random number
	void operator()( Arg &x)
	{
		x+=bound*( gsl_rng_uniform( gslrand1)-0.5);
	}
};

//*****************************************************************************
class Tsimmain
//-----------------------------------------------------------------------------
//! This is the main class of the Simulation encapsulating everything necessary to do
//! the job
//*****************************************************************************
{
	int sw; //!< which electron has been moved
	double bound;//!< the boundary for the random number distribution
	double alpha;//!< the parameter to be variated
	double q;//!< used for the EL calculation

	double r1, r2, r12; //! the distances
	double ro1, ro2, ro12; //!< the old distances
	double r112, r212;//!< scalar products of r1, r2 with r12
	double opsi2;//!< the old intensity
	double npsi2;//!< the new intensity

	double cEL;//!< cumulative EL value
	double csEL;//!< cumulative squared EL value

	vector<double>e1;//!< electron number one
	vector<double>e2;//!< electron number two
	vector<double>e12;//!< the distance vector
	vector<double>h;//!< if you know me you know what this does
	vector<double>old;//!< the old position of particle sw
	vector<double>conold;//!< the old connection vector
	
	vector<double>pcf;//!< the pair correlation function
	vector<double>::iterator vit;//!< vector iterators can allways come in handy

	gsl_rng *gslrand1; //!< the gsl random number generator 

	vAddRand<double> *ar;//!< instance of the unary_function to add a random number to a vector

	//! \brief performs one MC step
	//! \return if the step has been accepted
	bool MCstep();

	//! \brief calculate the distances
	inline void calcDistances();

	//! \brief calculate the local energy
	//! \return the local energy
	inline double EL();

	//! \brief calculate square of the wave function
	//! \return the square of the wave function
	inline double psi2();
public:

	//! \brief Constructor
	Tsimmain();

	//! \brief Destructor
	~Tsimmain();

	//! \brief function Initializing everythin necessary
	//! \param inN the dimension
	void Init();
	
	//! \brief this is the function doing the real work
	bool Execute();
	
	//! \brief this writes the pcf to pcf.dat
	void writePCF();
};

//*****************************************************************************
// Inline functions
//*****************************************************************************
//*****************************************************************************
void Tsimmain::calcDistances()
//*****************************************************************************
{
	std::transform( e1.begin(), e1.end(), e1.begin(), h.begin(), std::multiplies<double>());
	r1 = sqrt( std::accumulate (h.begin (), h.end (), 0.));
	
	std::transform( e2.begin(), e2.end(), e2.begin(), h.begin(), std::multiplies<double>());
	r2 = sqrt( std::accumulate (h.begin (), h.end (), 0.));
	
	std::transform( e2.begin(), e2.end(), e1.begin(), e12.begin(), std::minus<double>());
	std::transform( e12.begin(), e12.end(), e12.begin(), h.begin(), std::multiplies<double>());
	r12 = sqrt( std::accumulate (h.begin (), h.end (), 0.));
}

//*****************************************************************************
double Tsimmain::EL()
//*****************************************************************************
{
	q = 2.+2.*alpha*r12;
	std::transform(  e1.begin(), e1.end(), e12.begin(), h.begin(), std::multiplies<double>());
	r112 = std::accumulate (h.begin (), h.end (), 0.);
	std::transform(  e2.begin(), e2.end(), e12.begin(), h.begin(), std::multiplies<double>());
	r212 = std::accumulate (h.begin (), h.end (), 0.);

	return ( 1./(r12)-4.-4./( q*q*q*q)-4.*( r112/( r1*r12)-r212/( r2*r12))/( q*q)-8./( q*q*q*r12)); // keep an eye on that
}

//*****************************************************************************
double Tsimmain::psi2()
//*****************************************************************************
{
	double h4 = exp( ( r12/(2.+2.*alpha*r12))-2.*r1-2.*r2);
	return h4*h4;
}

#endif