/*************************************************************************************
* *
* ***** OpenThermo ***** *
* Calculation of thermodynamic functions from molecular data *
* Copyright 2008 Konstantin Tokarev <annulen@users.sourceforge.net> *
* and others *
* *
*************************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License. See COPYING for *
* more details *
* *
*************************************************************************************
* Module name : pes.hpp *
* Author : Tokarev. K *
* Last modification : 2008/07/23 *
* Description : This module contains declarations of class 'PES', *
* which is used for approximation and storage of potential *
* curves for internal rotors *
* *
*************************************************************************************/
#ifndef PES_HPP
#define PES_HPP
#include "thermo.h"
#include <iostream>
#include "hash.h"
#include <vector>
// PES - Potential Energy Surface
/// Type of PES approximation
enum PES_Approx_Type { FOURIER, ///< Fourier trigonometric polynom interpolation
POLYNOM, ///< Polynomial approximation
SPLINE ///< Spline interpolation
};
// Only 1D version now
/**
* \class PES pes.hpp
* \brief Representation of Potential Energy Surface
* \author Konstantin Tokarev
*
* The PES object is used for approximation and storage of potential surfaces
* for internal rotations. Main purpose of PES object is to return values
* of energy at arbitrary value of argument, encapsulating way this PES
* was defined (analytical or discreet) and interpolation technique in latter
* case. Currently only 1D PES (a.k.a. rotation barrier) is implemented
*/
class PES
{
public:
/**
* Constructor
* \param begin left border for coordinate change
* \param end right border for coordinate change
*/
PES (double begin, double end); // Borders for generic coordinate change
/// Destructor
~PES ();
// Access
double GetHeight () const; //!< Get barrier height (difference between max and min) in J/mol
double GetEnergy (double q) const; //!< Get energy at given coordinate \e q in J/mol
double dEnergy (double q) const; //!< Get energy 1st derivative at given coordinate \e q in J/mol
double d2Energy (double q) const; //!< Get energy 2nd derivative at given coordinate \e q in J/mol
const Hash<double,double> & GetPoints () const;
const std::vector<double> & GetMin () const;
const std::vector<double> & GetMax () const;
const std::vector<double> & GetCoeff () const;
const std::vector<double> & GetLevels () const; //!< Get energy levels \sa CalculateLevels
void SetSigma (int sigma); //!< Set rotation number of barrier
//bool GetSimple () const;
// Fast initialization
void SetPES (double h, double sigma);
void SetPES (Hash<double,double> & energy,
PES_Approx_Type a=FOURIER, int l=0);
/// Trigger calculation of energy levels
void CalculateLevels(const double I, const int N);
#ifdef OBJECT_COUNTERS
static int NumOfPES;
#endif
private:
void DoApprox ();
bool FourierApprox ();
//bool PolynomApprox ();
bool SplineApprox ();
void FindMinMax ();
void Newton (double &x);
double FindMin () const;
double FindMax () const;
void SetZeroX (double x);
void SetZeroY (double y);
//void Crowford (); // find levels
void CleanCoeff ();
Hash<double,double> itsEnergyPoints;
double itsBegin;
double itsEnd;
double itsHeight;
std::vector<double> itsLevel;
std::vector<double> itsFourierCoeff;
unsigned int itsL_max;
std::vector<double> itsPolynomCoeff;
std::vector<double> itsMax;
std::vector<double> itsMin;
PES_Approx_Type itsApprox;
bool isSimple;
};
typedef PES PES_1D;
#endif /* PES_HPP */