/*************************************************************************************
* *
* ***** 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 definitions of class 'PES', *
* which is used for approximation and storage of potential *
* curves for internal rotors *
* *
*************************************************************************************/
#include "pes.hpp"
#include "mathmodule.hpp"
const double Eps = 1e-10;
const double pi = 3.14159265358979323846;
PES::PES (double begin, double end)
{
itsApprox = FOURIER;
itsBegin = begin;
itsEnd = end;
itsL_max = 0;
}
PES::~PES ()
{
}
// Access
double PES::GetHeight ()
{
// if rotor is asymmetric, this value is meaningless!
// need for quantum solution
//if (itsHeight == 0) ...
return itsHeight;
}
double PES::GetEnergy (double q) const
{
double e = 0;
switch (itsApprox)
{
case FOURIER:
e = itsFourierCoeff.at(0)/2;
for (unsigned int l=1; l<=itsL_max; l++)
e += (itsFourierCoeff.at(2*l-1)*cos(l*q)+itsFourierCoeff.at(2*l)*sin(l*q));
break;
case POLYNOM:
for (int i=0; i<itsPolynomCoeff.size(); i++)
e += itsPolynomCoeff.at(i) * pow(q,i);
break;
case SPLINE:
break;
}
return e;
}
double PES::dEnergy (double q) const
{
double e = 0;
switch (itsApprox)
{
case FOURIER:
//MaxDegree
for (unsigned int l=1; l<=itsL_max; l++)
e += l*(itsFourierCoeff.at(2*l)*cos(l*q)-itsFourierCoeff.at(2*l-1)*sin(l*q));
break;
case POLYNOM:
for (int i=1; i<itsPolynomCoeff.size(); i++)
e += i*itsPolynomCoeff.at(i) * pow(q,i-1);
break;
case SPLINE:
break;
}
return e;
}
double PES::d2Energy (double q) const
{
double e = 0;
switch (itsApprox)
{
case FOURIER:
//MaxDegree
for (unsigned int l=1; l<=itsL_max; l++)
e -= l*l*(itsFourierCoeff.at(2*l)*sin(l*q)+itsFourierCoeff.at(2*l-1)*cos(l*q));
break;
case POLYNOM:
for (int i=1; i<itsPolynomCoeff.size(); i++)
e += i*(i-1)*itsPolynomCoeff.at(i) * pow(q,i-2);
break;
case SPLINE:
break;
}
return e;
}
const vector<double> & PES::GetMin () const
{
return itsMin;
}
const vector<double> & PES::GetMax () const
{
return itsMax;
}
const vector<double> & PES::GetCoeff () const
{
switch (itsApprox)
{
case FOURIER: return itsFourierCoeff; break;
case POLYNOM: return itsPolynomCoeff; break;
case SPLINE: break;
}
}
const map<double,double> & PES::GetPoints () const
{
return itsEnergyPoints;
}
void PES::SetPES (double h, double sigma)
{
itsEnergyPoints.clear();
itsLevels.clear();
CleanCoeff();
itsHeight = h;
itsApprox = FOURIER;
itsFourierCoeff.reserve(sigma+1);
itsFourierCoeff.push_back(h/2);
for (int i=0; i<sigma; i++)
itsFourierCoeff.push_back(0);
itsFourierCoeff.push_back(-h/2);
// Problem: alpha=0 for equil or it's real angle?
}
void PES::SetPES (map<double,double> & energy, PES_Approx_Type a, int l)
{
itsEnergyPoints.clear();
itsLevels.clear();
CleanCoeff();
itsEnergyPoints = energy;
itsApprox = a;
itsL_max = l;
DoApprox();
FindMinMax();
SetZeroY(GetEnergy(FindMin()));
CleanCoeff();
DoApprox();
}
void PES::CleanCoeff ()
{
switch (itsApprox)
{
case FOURIER: itsFourierCoeff.clear(); break;
case POLYNOM: itsPolynomCoeff.clear(); break;
case SPLINE: break;
}
}
void PES::DoApprox ()
{
switch (itsApprox)
{
case FOURIER:
if(FourierApprox())
break;
//case POLYNOM:
// if(PolynomApprox());
// break;
case SPLINE:
if(SplineApprox())
{
itsApprox = SPLINE;
break;
}
default:
throw GenericError("## ERROR : Can't approximate PES!");
}
}
bool PES::FourierApprox ()
{
double s=0;
int NNodes=0;
for (map<double,double>::const_iterator it=itsEnergyPoints.begin();
it!=itsEnergyPoints.end(); it++,NNodes++)
s += it->second; //y[k];
// NNodes MUST be even!!!
if (NNodes % 2 != 1)
return false;
if (itsL_max == 0)
itsL_max=(NNodes-1)/2;
//itsFourierCoeff // NNodes
itsFourierCoeff.push_back(2*s/NNodes);
for (int l=1; l<=itsL_max; l++)
{
double a = 0, b = 0;
//for (int k=0; k<NNodes; k++)
for (map<double,double>::const_iterator it=itsEnergyPoints.begin();
it!=itsEnergyPoints.end(); it++)
{
a += (it->second)*cos(l*(it->first));
b += (it->second)*sin(l*(it->first));
}
// a_l*cos(lx)
//itsFourierCoeff.at(2*l-1) += 2*a/NNodes;
itsFourierCoeff.push_back(2*a/NNodes);
// b_l*sin(lx)
itsFourierCoeff.push_back(2*b/NNodes);
}
return true;
}
bool PES::SplineApprox ()
{
// To be implemented
return false;
}
/*void PES::PrintPES (ostream & os, double step)
{
for (double q=itsBegin; q<itsEnd; q+=step)
os << setprecision(5) << setw(10) << q << "\t"
<< setprecision(12) << setw(14) << GetEnergy(q) << endl;
}*/
void PES::SetZeroX (double x)
{
map<double,double> newEnergyPoints;
for (map<double,double>::iterator it=itsEnergyPoints.begin(); it!=itsEnergyPoints.end(); it++)
newEnergyPoints[it->first-x] = it->second;
itsEnergyPoints = newEnergyPoints;
}
void PES::SetZeroY (double y)
{
for (map<double,double>::iterator it=itsEnergyPoints.begin(); it!=itsEnergyPoints.end(); it++)
it->second -= y;
}
void PES::FindMinMax ()
{
// In case of 1D PES minimums and maximums alternate
// That's why it's enough to find stationary points
// finding places where derivative is zero
// 1. Find "left" stationary poits from input
// 2. Find precise stationary points using Newton algorithm
//vector<double> minA, minB, maxA, maxB;
// Check first and last point
//int NNodes=0;
vector<double> x,y;
for (map<double,double>::const_iterator it=itsEnergyPoints.begin();
it!=itsEnergyPoints.end(); it++)
{
x.push_back(it->first);
y.push_back(it->second);
}
//for (unsigned int i=0; i<y.size(); i++)
// cerr << i<< "\t" << x.at(i) << "\t" << y.at(i) << endl;
/*int i=0;
map<double,double>::const_iterator it=itsEnergyPoints.begin();
it++;
i++;*/
//double min = it->second; //itsEnergyPoints[it->first];
/*itsMax.push_back(itsBegin);
itsMin.push_back(itsBegin);*/
//while (it!=itsEnergyPoints.end())
itsMin.clear();
itsMax.clear();
for (unsigned int i=0; i<y.size(); i++)
{
//double x = it->first;
//cout << x << endl;
/*if((itsEnergyPoints[i-1] < itsEnergyPoints[i])
&& (itsEnergyPoints[i] >= itsEnergyPoints[i+1]))
itsMax.push_back(x);
if((itsEnergyPoints[i-1] > itsEnergyPoints[i])
&& (itsEnergyPoints[i] <= itsEnergyPoints[i+1]))
itsMin.push_back(x);
it++;
i++;*/
if((y.at((i-1+y.size())%y.size()) <= y.at((i+y.size())%y.size())) && (y.at(i) > y.at((i+1)%y.size())))
itsMax.push_back(x.at(i));
if((y.at((i-1+y.size())%y.size()) >= y.at((i+y.size())%y.size())) && (y.at(i) < y.at((i+1)%y.size())))
itsMin.push_back(x.at(i));
}
//itsMax.push_back(itsEnd);
//itsMin.push_back(itsEnd);
for (unsigned int i=0; i<itsMin.size(); i++)
{
Newton(itsMin.at(i));
}
for (unsigned int i=0; i<itsMax.size(); i++)
{
Newton(itsMax.at(i));
}
/*
double x=itsBegin;
Newton(x);
for (int i=1; i<itsMin.size()-1; i++)
{
if (x == itsMin.at(i))
{
//delete 0 from Min and Max
goto BeginDeleted;
}
}
for (int i=1; i<itsMax.size()-1; i++)
{
if (x == itsMax.at(i))
{
//delete 0 from Min and Max
goto BeginDeleted;
}
}
BeginDeleted:
//*/
}
double PES::FindMin () const
{
double min = itsMin.at(0);
for (unsigned int i=1; i<itsMin.size(); i++)
{
if (GetEnergy(itsMin.at(i)) < GetEnergy(min))
min = itsMin.at(i);
}
return min;
}
void PES::Newton (double &x)
{
//cerr << x/pi*180 << "\t";
double dx=0;
do
{
//cerr << x << "\t" << dEnergy(x) << "\t" << d2Energy(x) << endl;
dx = dEnergy(x)/d2Energy(x);
x -= dx;
}while (fabs(GetEnergy(x)-GetEnergy(x+dx)) < Eps);
//cerr << x/pi*180 << endl;
}