#ifndef ROTOR_HPP
#define ROTOR_HPP
#include "thermo.h"
#include "pes.hpp"
#include "structure.hpp"
/**
* \enum ROTATION Type of rotation input data
*/
enum ROTATION {FREE, //!< Free rotation
HINDERED_HARMONIC, //!< Hindered rotation with barrier defined via Fourier coefficients
HINDERED_POLYNOM, //!< Hindered rotation with barrier defined via polynomical coefficients
HINDERED_DISCREET //!< Hindered rotation with barrier defined via xy plot
};
/**
* Calculates value of integral expression in partition function of internal
* rotor. Equation:\n
* \f$I=\int\limits_0^{2\pi}e^{-\frac {E(\phi)} {RT} }d\phi\f$
* \param T temperature in K
* \param Barrier PES object which provides \f$E(\phi)\f$
* \sa PES::GetEnergy
*/
double RotationIntegral(double T, PES_1D & Barrier);
/**
* \enum Rotation_Theory Theory used for calculation of partition function
* of rotor
*/
enum Rotation_Theory { CLASSIC, //!< Use classic equations
QUANTUM //!< Use energy levels of quantum rotor
};
/**
* \class Rotor rotor.hpp
* \brief Representation of internal rotor
* \author Konstantin Tokarev
*
* The Rotor class is a Structure subclass that provides an object for
* internal rotor, which may be a member of Molecule or another Rotor.
*
* \sa \ref rotors, \ref rotor
*/
class Rotor : public Structure
{
public:
/// Constructor
Rotor ();
/// Destructor
~Rotor ();
// Access methods
// Contents
void SetAtoms (int NAtoms);
void SetAtomNumbers (int * atoms);
int GetAtomNumber (int n) const;
// Connection with other part of molecule
/// Set atom from parent structure where rotor is connected
void SetSupport (Atom * theAtom);
/// Get atom from other part of parent structure where rotor is connected
Atom * GetSupport () const;
/**
* Set atom of rotor connected to other part of parent structure
* \param n number of head atom
*/
void SetHead (int n);
/**
* Set atom of rotor connected to other part of parent structure
* \param theAtom pointer to head atom
*/
void SetHead (Atom * theAtom);
/// Get atom of rotor connected to other part of parent structure
Atom * GetHead () const;
/// Set rotation number
void SetSigma (int N);
/// Set reduced inertial moment in 1st approximation \sa CalculateProjectionOfCenter
void SetI0 (double I0);
/// Get reduced inertial moment in 1st approximation \sa CalculateProjectionOfCenter
double GetI0 () const;
/// Set reduced inertial moment in 2nd approximation \sa CalculateInertialAxes
void SetIeff (double Ieff);
/// Get reduced inertial moment in 2nd approximation \sa CalculateInertialAxes
double GetIeff () const;
/// \return true if rotor is balanced
bool isBalanced() const { return itsBalanced; }
void SetRotType (ROTATION type);
ROTATION GetRotType () const;
Rotation_Theory GetTheory () const;
void SetTheory (Rotation_Theory theory);
void SetPointGroup (const char * pg);
/// Get the barrier of internal rotation
PES_1D & GetBarrier ();
/**
* Trigger calculation of coordinates of projection of mass center of rotor
* on its rotation axis
*/
void CalculateProjectionOfCenter ();
/**
* Trigger calculation of inertial axes of rotor, and define inernal
* coordinate system
* \warning Needs CalculateProjectionOfCenter to be called before
* \note There's no need to call this function for calculation of rediced moment
* if rotor is balanced
*/
void CalculateInertialAxes ();
// Stuff for calculation of reduced inertia momentum
double GetInn (Atom * A, Atom * B) const;
/**
* Get rotation moment of rotor about its rotation axis.
* Equation in coordinate system of rotor:\n
* \f$A = \sum\limits_i m_i(x_i^2+y_i^2)\f$
*/
double A() const; // moment about z
/**
* Get XZ product of rotor. Equation in coordinate system of rotor:\n
* \f$B = \sum\limits_i m_i x_i z_i\f$
* \sa CalculateInertialAxes
*/
double B() const; // xz product, sum xi*zi
/**
* Get YZ product of rotor. Equation in coordinate system of rotor:\n
* \f$C = \sum\limits_i m_i y_i z_i\f$
* \sa CalculateInertialAxes
*/
double C() const; // yz product, sum yi*zi
/**
* Get off-balance factor of rotor. Equation in coordinate system of rotor:\n
* \f$U = \sum\limits_i m_i x_i\f$
* \sa CalculateInertialAxes
*/
double U() const; // off-balance factor, sum mi*xi
/**
* Get the distance between the mass center of rotor and the mass center its parent.
* \sa GetMassCenter
*/
double r(unsigned int i) const;
/// Equivalent to CosineOfDirection(i, axis) \sa CosineOfDirection
double Alpha(int i, char axis) const { return CosineOfDirection(i, axis); }
double Beta(unsigned int i) const;
void backupCoordinates();
void Rotate(double alpha);
void restoreCoordinates();
/**
* Search for the point group of rotor. The only symmetry elements searched
* are axes passing through itsSupport
* \param LogName name of file to write log in
* \param exact if true, use strict criteria for search
*/
void TestSymmetry (const char * LogName, bool exact=false);
void LogAtoms (std::ostream & flog) const;
/// Print table with atoms to stream using coordinates from local coordinate system of rotor
void LogAtomsLocal (std::ostream & flog) const;
/// Print information about rotor: symmetry, moments, etc.
void LogInfo (std::ostream & flog) const;
private:
Atom * itsSupport;
Atom * itsHead;
int * itsAtomNumber;
double itsI0; //!< zeroth approximation of Ieff
double itsIeff;
double itsAlpha;
ROTATION itsRotType;
PES_1D itsBarrier;
Rotation_Theory itsRotationTheory;
bool itsBalanced;
double itsDeviation;
ColumnVector itsProjectionOfCenter;
void CalculateSigma ();
void GetProjectionOfCenter(double &x, double &y, double &z) const;
void GetLocalCoordinates(double &x, double &y, double &z) const;
};
#endif /* ROTOR_HPP */