#include "atom.hpp"
#include "rotor.hpp"
#include "control.hpp"
#include "exception.hpp"
#include "mathmodule.hpp"
#include "matrix.h"
#include "symmetry.hpp"
#include <cmath>
#include <iomanip>
#include <fstream>
#include <cstring>
using namespace std;
//extern const double Na, pi, c, h, R;
const double A2m = 1e-10; // Angstrom in meters
// Local declarations
double GenericRotor (double T, double alfa, PES_1D & Barrier);
// Local definitions
double GenericRotor (double T, double alfa, PES_1D & Barrier)
{
double E = Barrier.GetEnergy(alfa);
// energy per mol
return exp(-E/(R*T));
}
double RotationIntegral(double T, PES_1D & Barrier)
{
double I = Integrate(GenericRotor, 0, 2*pi, T, Barrier);
return I;
}
// Definitions of Rotor
Rotor::Rotor(): itsSupport(0), itsHead(0),
itsAtomNumber(0), itsI0(0), itsIeff(0), itsAlpha(0), itsRotType(HINDERED_POLYNOM),
itsBarrier(0,2*pi), itsRotationTheory(CLASSIC), itsBalanced(false),
itsProjectionOfCenter(3)
{
//cerr << "\nRotor constructor...";
}
Rotor::~Rotor()
{
//cerr << "\nRotor destructor...";
delete [] itsAtomNumber;
}
void Rotor::SetAtoms (int NAtoms)
{
itsAtomNumber = new int[NAtoms];
Structure::SetAtoms(NAtoms);
}
void Rotor::SetAtomNumbers (int * atoms)
{
for (int i=0; i<itsNAtoms; i++)
itsAtomNumber[i] = atoms[i];
}
int Rotor::GetAtomNumber (int n) const
{
return itsAtomNumber[n];
}
void Rotor::SetSupport (Atom * theAtom)
{
itsSupport = theAtom;
}
Atom * Rotor::GetSupport () const
{
return itsSupport;
}
void Rotor::SetHead (int n)
{
itsHead = itsAtom[n];
}
void Rotor::SetHead (Atom * theAtom)
{
itsHead = theAtom;
}
Atom * Rotor::GetHead () const
{
return itsHead;
}
void Rotor::SetSigma (int N)
{
if (N>1)
itsTopType = SYMMETRIC;
Structure::SetSigma(N);
itsBarrier.SetSigma(itsSigma);
}
void Rotor::SetI0 (double I0)
{
itsI0 = I0;
}
double Rotor::GetI0 () const
{
return itsI0;
}
void Rotor::SetIeff (double Ieff)
{
itsIeff = Ieff;
}
double Rotor::GetIeff () const
{
return itsIeff;
}
void Rotor::SetRotType (ROTATION type)
{
itsRotType = type;
}
ROTATION Rotor::GetRotType () const
{
return itsRotType;
}
Rotation_Theory Rotor::GetTheory () const
{
return itsRotationTheory;
}
void Rotor::SetTheory (Rotation_Theory theory)
{
itsRotationTheory = theory;
}
PES & Rotor::GetBarrier ()
{
return itsBarrier;
}
double Rotor::GetInn (Atom * A, Atom * B) const
{
// Returns moment for axis going through A and B
// This function doesn't use itsInertialTensor
double x1, y1, z1, x2, y2, z2, I;
A->GetCoordinates(x1, y1, z1);
B->GetCoordinates(x2, y2, z2);
I = 0;
for (int i=0; i<itsNAtoms; i++)
{
double xi, yi, zi, mi, hi2;
itsAtom[i]->GetCoordinates(xi,yi,zi);
mi = itsAtom[i]->GetWeightSI();
hi2 = (xi-x1)*(xi-x1) + (yi-y1)*(yi-y1) + (zi-z1)*(zi-z1) -
((x2-x1)*(xi-x1)+(y2-y1)*(yi-y1)+(z2-z1)*(zi-z1))*((x2-x1)*
(xi-x1)+(y2-y1)*(yi-y1)+(z2-z1)*(zi-z1)) /
((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
I += mi*hi2*A2m*A2m; // Convert A to m !!!
}
return I;
}
double Rotor::A() const
{
return GetInn(itsHead, itsSupport);
}
double Rotor::B() const
{
double b = 0;
double xi, yi, zi, mi;
for (int i=0; i<itsNAtoms; i++)
{
itsAtom[i]->GetCoordinates(xi,yi,zi);
GetLocalCoordinates(xi,yi,zi);
mi = itsAtom[i]->GetWeightSI();
b += mi*xi*zi*A2m*A2m;
}
return b;
}
double Rotor::C() const
{
double c = 0;
double xi, yi, zi, mi;
for (int i=0; i<itsNAtoms; i++)
{
itsAtom[i]->GetCoordinates(xi,yi,zi);
GetLocalCoordinates(xi,yi,zi);
mi = itsAtom[i]->GetWeightSI();
c += mi*yi*zi*A2m*A2m;
}
return c;
}
double Rotor::U() const
{
double u = 0;
double xi, yi, zi, mi;
for (int i=0; i<itsNAtoms; i++)
{
itsAtom[i]->GetCoordinates(xi,yi,zi);
GetLocalCoordinates(xi,yi,zi);
mi = itsAtom[i]->GetWeightSI();
u += mi*xi*A2m;
}
return u;
}
double Rotor::Beta(unsigned int i) const
{
return Alpha(i, 'z')*A() - Alpha(i, 'x')*B() - Alpha(i, 'y')*C()
+ U()*(Alpha(indexRing(i-1, 3), 'y')*r(indexRing(i+1, 3))
- Alpha(indexRing(i+1, 3), 'y')*r(indexRing(i-1, 3)));
}
double Rotor::r(unsigned int i) const
{
if (!itsParent)
throw GenericError("Rotor::r requires parent!");
double Xc, Yc, Zc, Xp, Yp, Zp;
GetMassCenter(Xc, Yc, Zc);
itsParent->GetMassCenter(Xp, Yp, Zp);
//cerr << "parent mass ceter: " << Xp << " " << Yp << " " << Zp << endl;
switch(i)
{
case 0: return (Xc-Xp)*A2m;
case 1: return (Yc-Yp)*A2m;
case 2: return (Zc-Zp)*A2m;
default:
cerr << "i=" << i << endl;
throw GenericError("Rotor::r - incorrect coordinate given!");
}
}
void Rotor::Rotate(double alpha)
{
// Rotate around axis of Rotor
double x1, y1, z1, x2, y2, z2, x, y, z, l;
long double alpha_change = alpha - itsAlpha;
Matrix R(3,3);
ColumnVector V1(3);
ColumnVector V2(3);
itsSupport->GetCoordinates(x1, y1, z1);
itsHead->GetCoordinates(x2, y2, z2);
l = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
x = (x2-x1)/l;
y = (y2-y1)/l;
z = (z2-z1)/l;
// Matrix of rotation
R(1,1) = cos(alpha_change)+(1-cos(alpha_change))*x*x;
R(1,2) = (1-cos(alpha_change))*x*y-sin(alpha_change)*z;
R(1,3) = (1-cos(alpha_change))*x*z+sin(alpha_change)*y;
R(2,1) = (1-cos(alpha_change))*y*x+sin(alpha_change)*z;
R(2,2) = cos(alpha_change)+(1-cos(alpha_change))*y*y;
R(2,3) = (1-cos(alpha_change))*y*z-sin(alpha_change)*x;
R(3,1) = (1-cos(alpha_change))*z*x-sin(alpha_change)*y;
R(3,2) = (1-cos(alpha_change))*z*y+sin(alpha_change)*x;
R(3,3) = cos(alpha_change)+(1-cos(alpha_change))*z*z;
for (int i=0; i<itsNAtoms; i++)
{
double xi, yi, zi;
itsAtom[i]->GetCoordinates(xi,yi,zi);
V1 << xi << yi << zi;
// Transformation of coordinates
V2 = R*V1;
itsAtom[i]->SetCoordinates(V2(1),V2(2),V2(3));
}
itsAlpha = alpha;
}
void Rotor::SetPointGroup (const char * pg)
{
Structure::SetPointGroup(pg);
CalculateSigma();
}
void Rotor::TestSymmetry (const char * LogName, bool exact)
{
char pg[15], pg_brutal[15];
double support[3];
char C1[] = "C1";
//#ifdef POSIX
string tmpname = tmpnam(0);
//#else
// char tmpname[] = tmpnam(0);
//#endif
SymmetryState s;
ofstream tmp(tmpname.c_str());
ofstream fout;
fout.open(LogName, ios::app);
// Preparing input file for Symmetry program
ExportCoordinates(tmp);
tmp.close();
/*// Find number of Head
for (int i=0; i<itsNAtoms; i++)
if (itsAtom[i] == itsHead)
{
head = i;
break;
}
// -1 - exception*/
// Running Symmetry
fout << "\n***** Symmetry (Brute force symmetry analyzer) Author: S. Patchkovskii, 1996,2003 *****" << endl;
fout << "\nStarting Symmetry for rotor...\n" << endl;
fout << "Search thresholds\n\tPrimary: " << DefaultSymmetryPrimary
<< "\n\tFinal: " << DefaultSymmetryFinal
<< "\n\tMaxit: " << DefaultSymmetryMaxit << endl;
fout.close();
itsSupport->GetCoordinates(support[0], support[1], support[2]);
s = symmetry_of_axis(pg, support, tmpname.c_str(), LogName);
fout.open(LogName, ios::app);
switch (s)
{
case OK:
if (exact != true)
{
fout << "\nRestarting Symmetry to find higher axis...\n" << endl;
fout << "Search thresholds\n\tPrimary: " << Control::instance()->SymmetryPrimary()
<< "\n\tFinal: " << Control::instance()->SymmetryFinal()
<< "\n\tMaxit: " << Control::instance()->SymmetryMaxit() << endl;
fout.close();
s = symmetry_of_axis(pg_brutal, support, tmpname.c_str(), LogName, -2,
Control::instance()->SymmetryPrimary(),
Control::instance()->SymmetryFinal(),
Control::instance()->SymmetryMaxit());
fout.open(LogName, ios::app);
switch (s)
{
case OK:
SetPointGroup(pg_brutal);
fout << "\nSetting " << pg_brutal << " point group" << endl;
if(strcmp(pg, pg_brutal))
{
fout << "WARNING! " << pg_brutal << " isn't exact point group!" << endl; //\n You can choose
cerr << "# WARNING! Setting " << pg_brutal << ", but it isn't exact point group of rotor" << endl; //\n You can choose
}
goto cleantmp;
case UNKNOWN_GROUP:
fout << "\nSymmetry program can't determine point group, but rotor might have axis of higher order" << endl;
// falling down
default:
SetPointGroup(pg);
fout << "\nSetting " << pg << " point group" << endl;
goto cleantmp;
}
}
else
{
SetPointGroup(pg);
fout << "\nSetting " << pg << " point group" << endl;
goto cleantmp;
}
case UNKNOWN_GROUP:
// exception!
fout << "\nSymmetry program can't determine point group, but rotor might have axis of higher order" << endl;
goto cleantmp;
default:
fout << "Can't recognize point group, C1 is used as default" << endl; //\nUse keyword \'sym\' to specify point group" << endl;
cerr << "# WARNING! Can't recognize point group, C1 is used default" << endl;
SetPointGroup(C1);
goto cleantmp;
}
cleantmp:
fout << "**********" << endl;
tmp.open(tmpname.c_str());
tmp.close();
}
void Rotor::LogAtoms (ostream & flog) const
{
flog << "\n\t\t********** Rotor XYZ coordinates **********\n";
Structure::LogAtoms(flog);
flog << "Head atom:\n" << *itsHead;
flog << "Support atom:\n" << *itsSupport;
flog << endl;
}
void Rotor::LogAtomsLocal (ostream & flog) const
{
flog << itsNAtoms << endl;
//flog << "Atom\t X\t Y\t Z\t M\t Spin" << endl;
double xi,yi,zi;
for (int i=0, k=1; i<itsNAtoms; i++, k++) {
itsAtom[i]->GetCoordinates(xi,yi,zi);
GetLocalCoordinates(xi,yi,zi);
flog << itsAtom[i]->GetName() << " " << xi << " " << yi
<< " " << zi << endl;
}
flog << endl;
}
void Rotor::LogInfo (ostream & flog) const
{
vector<double> coeff, min, max;
flog << "\n***** Rotor Summary *****" << endl;
flog << "Point group: " << itsPointGroup << endl;
switch (itsTopType)
{
case LINEAR:
flog << "Linear top";
break;
case SYMMETRIC:
flog << "Symmetric top";
break;
case ASYMMETRIC:
flog << "Asymmetric top";
break;
default:
flog << "Unknown top\n ##ERROR! Something went wrong!";
break;
}
if (itsBalanced)
{
flog << "\nBalanced top";
flog << "\n Deviation of mass center from rotation axis: " << itsDeviation << " ang";
}
else
{
flog << "\nNon-babalnced top";
flog << "\n Deviation of mass center from rotation axis: "
<< itsDeviation << " ang";
flog << "\n Off-balance factor: "
<< U()*1e5 << " [g*cm]\t" << U()/(1e-3/Na)/A2m << " [a.u.*ang]";
}
flog << "\n I(rotor) = " << A()*1e7 << " [g*cm^2]";
flog << "\n Ired = " << itsI0*1e7 << " [g*cm^2] (1st approximation)";
flog << "\n Ired = " << itsIeff*1e7 << " [g*cm^2] (2nd approximation)";
flog << "\nRotation number: " << itsSigma << endl;
switch (itsRotType)
{
case FREE:
flog << "Free rotation" << endl;
break;
case HINDERED_POLYNOM:
break;
case HINDERED_HARMONIC:
case HINDERED_DISCREET:
flog << "Hindered rotation" << endl;
flog << "\nFourier coefficients:\nl\t"<< setw(12) << setiosflags(ios::right)
<< "cos" << "\t" << setw(12) << setiosflags(ios::right) << "sin\n";
coeff = itsBarrier.GetCoeff();
flog << 0 << "\t" << setw(12) << setiosflags(ios::showpoint) << setiosflags(ios::right) << coeff.at(0) << endl;
for (unsigned i=1; i<=(coeff.size()-1)/2; i++)
{
flog << i << "\t" << setw(12) << setiosflags(ios::showpoint) << setiosflags(ios::right) << coeff.at(2*i-1)
<< "\t" << setw(12) << setiosflags(ios::showpoint) << setiosflags(ios::right) << coeff.at(2*i) << endl;
}
flog << endl;
/*flog << "\nApproximation results:\n";
for (int x=0; x<360; x+=15)
{
flog << "\t" << x << "\t";
flog.precision(12);
flog << itsBarrier.GetEnergy(x*pi/180)/1000 << endl;
}*/
if (itsRotationTheory == QUANTUM)
{
flog << "Energy levels\n \t" << setw(12) << setiosflags(ios::right)
<< "[kJ/mol]\t" << setw(12) << setiosflags(ios::right) << "[1/cm]\n";
vector<double> l = itsBarrier.GetLevels();
unsigned int i=0;
do
{
flog << i << "\t" << setw(12) << setiosflags(ios::showpoint)
<< setiosflags(ios::right) << l.at(i)*Na/1000 << "\t" << setw(12)
<< setiosflags(ios::showpoint) << setiosflags(ios::right)
<< l.at(i)/(h*c) << "\n";
i++;
if ((i==21) && (i<l.size()-11))
{
i=l.size()-10;
flog << "...\t...\n";
}
}while (/*(l.at(i) <= itsBarrier.GetHeight()) &&*/ (i<l.size()-1));
flog << endl;
flog << "Overtone frequencies [1/cm]:\n";
for (int i=1; i<20; i++)
flog << i-1 << "->" << i << "\t" << (l.at(i)-l.at(i-1))/(h*c) << endl;
}
if (itsRotType == HINDERED_DISCREET)
{
flog << "\nApproximation results:\n";
for (int x=0; x<360; x+=15)
{
flog << "\t" << x << "\t";
flog.precision(12);
flog << itsBarrier.GetEnergy(x*pi/180)*1e-3 << endl;
}
/*min = itsBarrier.GetMin();
flog << "Harmonic frequencies in minima:\n";
for (unsigned i=0; i<min.size(); i++)
flog << "Angle: " << min.at(i)/pi*180 << "\tVibr. frequency: "
<< sqrt(itsBarrier.d2Energy(min.at(i))/(Na*itsIeff))/(2*pi*c) << " cm^-1" << endl;*/
min = itsBarrier.GetMin();
max = itsBarrier.GetMax();
flog << "\n\tx\tEnergy[kJ/mol]";
flog << "\nMinima:\n";
for (unsigned int i=0; i<min.size(); i++)
flog << i+1 << "\t" << min.at(i)/pi*180 << "\t" << itsBarrier.GetEnergy(min.at(i))*1e-3 << endl;
flog << "\nMaxima:\n";
for (unsigned int i=0; i<min.size(); i++)
flog << i+1 << "\t" << max.at(i)/pi*180 << "\t" << itsBarrier.GetEnergy(max.at(i))*1e-3 << endl;
flog << "\nMaximum Height[kJ/mol]:\t" << itsBarrier.GetHeight()*1e-3 << endl;
}
break;
default:
flog << "Unknown rotation\n ##ERROR! Something went wrong!" << endl;
}
flog << endl;
}
void Rotor::CalculateSigma ()
{
Structure::CalculateSigma();
if (itsSigma == 1)
itsTopType = ASYMMETRIC;
else
itsTopType = SYMMETRIC;
itsBarrier.SetSigma(itsSigma);
}
void Rotor::GetProjectionOfCenter(double &x, double &y, double &z) const
{
x = itsProjectionOfCenter(1);
y = itsProjectionOfCenter(2);
z = itsProjectionOfCenter(3);
}
void Rotor::CalculateProjectionOfCenter ()
{
ColumnVector A,B,C,I1(3),I2(3),I3(3);
itsSupport->GetCoordinates(A);
itsHead->GetCoordinates(B);
GetMassCenter(C);
I3 = B - A;//itsInertialAxes.Column(1);
I3 = I3/I3.NormFrobenius();
itsProjectionOfCenter = I3*DotProduct(I3,C-A) + A;
itsInertialAxes.Column(3) = I3;
itsDeviation = (itsProjectionOfCenter-C).NormFrobenius();
//cerr << "d = " << (itsProjectionOfCenter-C).NormFrobenius() << endl;
if(itsDeviation <= Control::instance()->BalancedRotor())
itsBalanced = true;
else
itsBalanced = false;
}
void Rotor::CalculateInertialAxes ()
{
// Requires CalculateProjectionOfCenter to be already called
ColumnVector A,B,C,I1(3),I2(3);
ColumnVector I3 = itsInertialAxes.Column(3);
GetMassCenter(C);
I1 = C - itsProjectionOfCenter; //itsInertialAxes.Column(3);
I1 = I1/I1.NormFrobenius();
// y is z*x
I2(1) = -I1(2)*I3(3) + I1(3)*I3(2);
I2(2) = -I1(3)*I3(1) + I1(1)*I3(3);
I2(3) = -I1(1)*I3(2) + I1(2)*I3(1);
I2 = I2/I2.NormFrobenius(); // prevent weird cases
itsInertialAxes.Column(1) = I1; //I1/I1.NormFrobenius();
itsInertialAxes.Column(2) = I2; //I2/I2.NormFrobenius();
//itsInertialAxes.Column(3) = I3; //I3/I3.NormFrobenius();
//itsInertialAxes.Column(1) = (C - itsProjectionOfCenter);
//itsInertialAxes.Column(1)
#ifdef DEBUG
cerr << "itsProjectionOfCenter:\n" << itsProjectionOfCenter << endl;
cerr << "Mass Center:\n" << C << endl;
cerr << "itsInertialAxes:\n" << itsInertialAxes << endl;
cerr << "itsInertialAxes.Determinant=" << itsInertialAxes.Determinant() << endl;
//if (itsInertialAxes.Determinant() != 1.0) {
// throw (GenericError("itsInertialAxes.Determinant != 1"));
// }
#endif
}
void Rotor::GetLocalCoordinates(double &x, double &y, double &z) const
{
ColumnVector e(3);
e << x << y << z;
// Translation
e -= itsProjectionOfCenter;
// Rotation
const ColumnVector internal = itsInertialAxes.t()*e;
// Return local coordinates
x = internal(1);
y = internal(2);
z = internal(3);
}