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

//! \file tsimbox.hpp
//! \date 15-11-2004

//! \class TParticle
//! \brief The Container class for one particle

//! \class TSimBox
//! \brief this class represents one simulation box 

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

//     This file is part of gibbsens.
//
//     gibbsens 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.
//
//     gibbsens 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 gibbsens; if not, write to the Free Software
//     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

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

#ifndef TSIMBOX_HPP
#define TSIMBOX_HPP

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
// #include <gsl/gsl_fft_real.h>
// #include <gsl/gsl_fft_complex.h>
// #include <gsl/gsl_errno.h>

#include <list>
#include <string>
#include <vector>
#include <algorithm>
#include <numeric>
#include <cmath>
#include <ctime>
#include <iostream>
#include <fstream>
#include <sstream>


const int adjsteps=1000;//!< number of steps before adjustment in the initialisation phase
// const double Vadjsteps=100;//!< number of Vsteps before adjustment
const bool debug=false;//!< debug output?
// const double meltmax=100.;//!< at which melting factor value we consider the stuff molten
const double startstep=0.9;//!< trial stepwith to start with
const double startvstep=0.01;//!< trial volume stepwith to start with
const double chstep=0.01;//!< rate at which the trial stepwitch will be changed
const double chvstep=0.8;//!< rate at which the trial volume stepwitch will be changed

const double rkb=3.16669e-6;//!< Boltzmann's constant in AU

class TParticle
{
	double a, b, c;
public:
	double pos[3];//!< the position of the particle

	//! \brief constructor
	//! \param x the x coordinate
	//! \param y the y coordinate
	//! \param z the z coordinate
	TParticle( const double x=0., const double y=0., const double z=0.){ pos[0]=x; pos[1]=y; pos[2]=z;}

	//! \brief copy constructor
	//! \param p the particle to copy
	TParticle( const TParticle *p){ pos[0]=p->pos[0]; pos[1]=p->pos[1]; pos[2]=p->pos[2];}


	//! \brief calculate the distance betweent this particle and a given one
	//! \param g The given particle
	//! \param L the current Box Length
	//! \return the distance
	double dist( const TParticle *g, const double L)
	{
		a=fabs( g->pos[0]-pos[0]);
		b=fabs( g->pos[1]-pos[1]);
		c=fabs( g->pos[2]-pos[2]);
		if( a>L/2.)a-=L;
		if( b>L/2.)b-=L;
		if( c>L/2.)c-=L;
		return sqrt( a*a+b*b+c*c);
	}

	//! \brief add something to the current position
	//! \param x the x coordinate
	//! \param y the y coordinate
	//! \param z the z coordinate
	//! \param L the current Box Length
	void addpos( const double x, const double y, const double z, const double L)
	{ 
		pos[0]=fmod( pos[0]+x+2.*L, L); 
		pos[1]=fmod( pos[1]+y+2.*L, L); 
		pos[2]=fmod( pos[2]+z+2.*L, L);
	}

	//! \brief scale all coordinates by a factor
	//! \param f the factor
	void scale( const double f)
	{
		pos[0]*=f;
		pos[1]*=f;
		pos[2]*=f;
	}
};

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


class TSimBox
{

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

	TSimBox *oBox; //!< Pointer to the other Box

	std::string name;//!< Name of this Box
	bool calcgofr;//!< wether or not to calculate the Pair Correlation Function
	double cells;//!< number of Unit Cells in one direction in the Initial Config
	double boxlen;//!< the Box Length (can change due to the Volume change)
	double halfbox;//!< half box length
	double m;//!< nhe Mass of the particles
	int N;//!< nhe Number of Particles currently in the Box (can change too)
	double T;//!< the Temperature in units of kB
	double beta;//!< the inverse Temperature in units of 1/kB
	double meltmax;//!< at which melting factor value we consider the stuff molten

//! Analysation Parameters
//! PCF stuff
	double pcfupper;//!< upper boundary for the pair correlation function
	int pcfrez;//!< number of bins for the pair correlation function
	double pcfdr;//!< binsize for the pair correlation function
	bool pcflog;//!< log scale for the pair correlation function
	std::vector<double> pcf;//!< the PCF Histogram
	std::vector<double>::iterator pcfit;//!< the PCF Histogram iterator

// Monte Carlo Parameters
	int ap;//!< the particle to relocate in the mc step
	double d;//!< the maximum displacement
	double dx;//!< the x displacement
	double dy;//!< the y displacement
	double dz;//!< the z displacement
	double E;//!< the total potential energy
	double Et;//!< the total potential energy of the trial step
	double br;//!< the distance buffer
	long mca;//!< the mc acceptance
	long steps;//!< the number of steps done in total
	long mcsteps;//!< the number of mc steps done in total
	double mcar;//!< the mc acceptance ratio

	TParticle buff;//!< buffer for the mc steps

//! volume step related parameters
	double V1; //!< current volume this box
	double V2; //!< current volume other box
	double V; //!< total volume
	double lnV; //!< log volume change
	double nV1; //!< new volume this box
	double nV2; //!< new volume other box
	double vd; //!< the volume displacement
	double VEn1; //!< the volume energy of box 1
	double VEn2; //!< the volume energy of box 2
	double f1; //!< the volume change factos
	double arg; // the acceptance exponent
	long va;//!< the v step acceptance
	double var;//!< the v step acceptance ratio
	long vsteps;//!< the number of v steps done in total

//! swapping step related stuff
	int which; //!< which particle to swap
	double NEn1; //!< the particle swap energy of box 1
	double NEn2; //!< the particle swap energy of box 2
	long sa;//!< the swap step acceptance
	double sar;//!< the swap step acceptance ratio
	long ssteps;//!< the number of swap steps done in total

// Silvera Goldstein Potential related stuff
// parameters
	double palpha; //! Silvera Golstein Potential parameter
	double pbeta; //! Silvera Golstein Potential parameter
	double pgamma; //! Silvera Golstein Potential parameter
	double pc6; //! Silvera Golstein Potential parameter
	double pc8; //! Silvera Golstein Potential parameter
	double pc9; //! Silvera Golstein Potential parameter
	double pc10; //! Silvera Golstein Potential parameter
	double prc; //! Silvera Golstein Potential parameter
	double pgammap; //! Silvera Golstein Potential parameter
	double pc6p; //! Silvera Golstein Potential parameter
	double pc8p; //! Silvera Golstein Potential parameter
	double pc9p; //! Silvera Golstein Potential parameter
	double pc10p; //! Silvera Golstein Potential parameter

// variables
	double pe; //! Silvera Golstein Potential variable
// 	double rdedr;
	double onr, onr2, onr6; //! Silvera Golstein Potential variable
	double px, py; //! Silvera Golstein Potential variable
	double pf; //! Silvera Golstein Potential variable
// 	double rdfdr;
	double pg; //! Silvera Golstein Potential variable
// 	double rdgdr;


//! vectors and matrices
	std::vector<TParticle> part;//!< List containing all the particles in the box
// 	std::vector<TParticle> buff;//!< buffer for the mc steps
	std::vector<TParticle>::iterator piti;//!< the particle iterator
	std::vector<TParticle>::iterator pitj;//!< the particle iterator

	std::list<std::list<double> > Em; //!< the energy matrix
	std::list<double>::iterator emitj; //!< the energy matrix j iterator
	std::list<std::list<double> >::iterator emiti; //!< the energy matrix i iterator
	std::list<double> Ebuff; //!< the Energy matrix line buffer
	std::list<double>::iterator Ebit; //!< the Energy matrix line buffer iterator

	std::list<std::list<double> > dist; //!< the distance matrix
	std::list<double>::iterator ditj; //!< the distance matrix j iterator
	std::list<std::list<double> >::iterator diti; //!< the distance matrix i iterator
	std::list<double> dbuff; //!< the distance matrix line buffer
	std::list<double>::iterator dbit; //!< the distance matrix line buffer iterator

	std::ostream *posout;//!< position output stream
	std::ostream *statout;//!< stats output stream

	//! \brief make the start lattice
	void makeStartCfg();

	//! \brief load the start lattice
	//! \param fname the name of the file to load
	//! \return true if everything worked
	bool loadCfg( std::string fname="");

	//! \brief calculate the energy matrix
	//! \param which which particle needs to be recalculated (-1 if the whole matrix gets redone)
	void calcEMatrix( const int which = -1);

	//! \brief update the energy matrix
	//! \param which which particle needs to be updated
	void updateEMatrix( const int which);


	//! \brief calculate the potential with a given |r|
	//! \param r the distance between the two particles
	//! \return the energy
	inline double EPot( const double r);

	//! \brief adjust the current displacement step
	//! \return is the step allready good?
	inline bool stepAdjust();

protected:
//! particle step related stuff

	//! \brief receive a particle
	//! \param p the particle
	//! \return the energy of the new config
	inline double recPart( const TParticle *p);

	//! \brief receive a particle (when box is empty)
	//! \return the energy of the new config
	inline double recPart();

	//! \brief accept the Particle swap step
	//! \param En the energy of the new configuration
	inline void accPart( const double En);

//! volume step related stuff

	//! \brief receive Volume
	//! \param 	f the volume to add
	//! \return the energy of the changed configuration
	inline double recVol( const double f);

	//! \brief accept the Volume change step
	//! \param En the energy of the new configuration
	inline void accVol( const double En);

	//! \brief calculate the energy of the new Volume configuration
	//! \param f the multiplicative factor
	//! \return the energy
	inline double getVolEn( const double f);
public:

	//! \brief standard constructor
	//! \param n name of this box ( used mainly for output)
	//! \param mass mass of the particles ( standard = 1)
	//! \param seed random number generator seed
	TSimBox( std::string n, const double mass =1., const int seed=0);
	//! \brief standard destructor standard 
	~TSimBox();

	//! \brief initialize the Box
	//! \param Box pointer to the other Box
	//! \param Num the number of particles to start with (needs to be 4*n^3 with n being int)
	//! \param bl the boxlength to start with
	//! \param Temp the Temperature in units of kB
	void Init( TSimBox *Box, const int Num, const double bl, const double Temp);

	//! \brief initialize the Box from a File
	//! \param Box pointer to the other Box
	//! \param fname the name of the file to read
	void Init( TSimBox *Box, std::string fname="");

	//! \brief setup the PCF calculation
	//! \param upper upper boundary of the PCF
	//! \param rez resolution of the PCF ( number of bins)
	//! \param logdist logarithmically equidistant ( not used)
	//! \todo write logscale stuff
	void initPCF( const int rez, const double upper=-1., const bool logdist=false);

	//! \brief Finalize the PCF calculation
	//! \param ncfg number of configs used for the PCF (std = internal steps counter)
	//! \param fname the filename to write the PCF to
	void finalizePCF( std::string fname = "PCF.dat", int ncfg = -1);

	//! \brief write the Structure Factor
	//! \param fname filename to write to
	//! \attention use after finalizePCF() ONLY!
	void writeSFactor( std::string fname = "Sofk.dat");

	//! \brief add the current step to the average
	inline void calcPCF();

	//! \brief calculate the Melt Factor
	//! \return the melt factor
	double getMelt();

	//! \brief do a Monte Carlo step
	//! \return acceptance
	bool mcStep();

	//! \brief do a particle swap step
	//! \return acceptance
	bool partStep();

	//! \brief do a volume step
	//! \return acceptance
	bool volStep();

	//! \brief save all important parameters
	//! \param fname the name of the file to save to
	void saveCfg( std::string fname="");

	//! \brief sets the position output stream
	inline void setPosOut( std::ostream *out = &std::cout);

	//! \brief sets the stats output stream
	inline void setStatOut( std::ostream *out = &std::cout);

	//! \brief prints out the current positions
	inline void printPos( const bool xmmol=false);

	//! \brief prints out the current Stats
	inline void printStats();

	//! \brief returns the volume
	//! \return the volume
	inline double getVol() const;

	//! \brief returns the Energy
	//! \return the Energy
	inline double getEn() const;

	//! \brief returns the Particle Number
	//! \return the Particle Number
	inline double getN() const;

	//! \brief adjust all steps to regulate the acceptance probability
	inline void adjustAllSteps();

// 	void printEm()
// 	{
// 		std::cout<<std::endl<<std::endl;
// 		for( emiti= Em.begin(); emiti!=Em.end(); ++emiti)
// 		{
// 			for( emitj=emiti->begin(); emitj!=emiti->end(); ++emitj)std::cout<<*emitj<<"  ";
// 			std::cout<<std::endl;
// 		}
// 	}

	void printDist()
	{
		std::cout<<std::endl<<std::endl;
		for( diti= dist.begin(); diti!=dist.end(); ++diti)
		{
			for( ditj=diti->begin(); ditj!=diti->end(); ++ditj)std::cout<<*ditj<<"  ";
			std::cout<<std::endl;
		}
	}
};

//*****************************************************************************
// inline implementations
//*****************************************************************************

//*****************************************************************************
void TSimBox::setPosOut( std::ostream *out)
//*****************************************************************************
{
	posout=out;
	*posout<<N<<std::endl;
}

//*****************************************************************************
void TSimBox::setStatOut( std::ostream *out)
//*****************************************************************************
{
	statout=out;
	*statout<<"Name of the Box =\t"<<name<<std::endl;
	*statout<<"Temperature =\t"<<T<<std::endl;
	*statout<<"Particle Mass =\t"<<m<<std::endl;
	*statout<<"Adjustment Steps =\t"<<adjsteps<<std::endl<<std::endl;

	*statout<<"Nonconstant parameters"<<std::endl;
	*statout<<"Steps\tE\tV\tN\trho\tAcc\tDis\tVAcc\tVDis\tPAcc\tfmelt"<<std::endl;
}

//*****************************************************************************
void TSimBox::printPos( const bool xmmol)
//*****************************************************************************
{
// 	*posout<<N<<std::endl;
	for( piti = part.begin(); piti!=part.end(); ++piti)
	{
		if( xmmol)*posout<<"H\t";
		if( !xmmol)*posout<<piti->pos[0]<<"\t"<<piti->pos[1]<<"\t"<<piti->pos[2];
		else *posout<<piti->pos[0]/boxlen<<"\t"<<piti->pos[1]/boxlen<<"\t"<<piti->pos[2]/boxlen;
		if( xmmol)*posout<<"\t0\t0\t0";
		*posout<<std::endl;
	}
}

//*****************************************************************************
void TSimBox::printStats()
//*****************************************************************************
{
	*statout<<steps<<" \t "<<getEn()<<" \t "<<getVol()<<" \t "<<getN()<<"\t"<<getN()/getVol()<<" \t "<<mcar<<" \t "<<d<<"\t"<<var<<"\t"<<vd<<" \t "<<sar<<" \t "<<getMelt()<<std::endl;
}

//***************************************************************************
void TSimBox::calcPCF()
//***************************************************************************
{
	int i,j, chn;
	if( N>=2)
	{
		for( diti=dist.begin(), i=0; i<N; ++diti, ++i)
			for( ditj=diti->begin(), j=i+2; j<N; ++ditj, ++j)
			{
				chn=static_cast<int>( *ditj/pcfdr);
				if( ( chn<pcfrez) && ( *ditj!=0.))pcf[chn]+=2.;
			}
		}
}

//***************************************************************************
double TSimBox::getVol() const
//***************************************************************************
{
	return boxlen*boxlen*boxlen;
}

//***************************************************************************
double TSimBox::getEn() const
//***************************************************************************
{
	return E;
}

//***************************************************************************
double TSimBox::getN() const
//***************************************************************************
{
	return N;
// 	return Em.begin()->size();
}

//***************************************************************************
void TSimBox::adjustAllSteps()
//***************************************************************************
{
//! volume step adjustment
	var=static_cast<double>( va)/static_cast<double>( vsteps);
// 	std::cerr<<va<<"\t"<<vsteps<<"\t"<<steps<<"\t"<<name<<std::endl;
 	if( var>0.55)
	{
		vd/=chvstep;
	}else if( var<0.45)
	{
		vd*=chvstep;
	}
	vsteps=0;
	va=0;

	//! mc step adjustment
	mcar=static_cast<double>( mca)/static_cast<double>( mcsteps);
	if( mcar>0.55)
	{
		d+=chstep;
	}else if( mcar<0.45)
	{
		d-=chstep;
	}
	mcsteps=0;
	mca=0;

	sar=static_cast<double>( sa)/static_cast<double>( ssteps);
	sa=ssteps=0;
}
#endif //TSIMBOX_HPP