source: src/molecule.hpp@ ecb799

/** \file molecule.hpp
 *
 * Class definitions of atom and molecule, element and periodentafel
 */

#ifndef MOLECULES_HPP_
#define MOLECULES_HPP_

using namespace std;

/*********************************************** includes ***********************************/

// GSL headers
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_heapsort.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_randist.h>

//// STL headers
#include <map>
#include <set>
#include <deque>
#include <list>
#include <vector>

#include <string>

#include "defs.hpp"
#include "graph.hpp"
#include "stackclass.hpp"
#include "tesselation.hpp"
#include "Patterns/Observer.hpp"
#include "Patterns/Cacheable.hpp"

/****************************************** forward declarations *****************************/

class atom;
class bond;
class BondedParticle;
class BondGraph;
class element;
class ForceMatrix;
class LinkedCell;
class molecule;
class MoleculeLeafClass;
class MoleculeListClass;
class periodentafel;
class Vector;

/******************************** Some definitions for easier reading **********************************/

#define MoleculeList list <molecule *>
#define MoleculeListTest pair <MoleculeList::iterator, bool>

#define DistancePair pair < double, atom* >
#define DistanceMap multimap < double, atom* >
#define DistanceTestPair pair < DistanceMap::iterator, bool>


/************************************* Class definitions ****************************************/

/** Structure to contain parameters needed for evaluation of constraint potential.
 */
struct EvaluatePotential
{
  int startstep;              //!< start configuration (MDStep in atom::trajectory)
  int endstep;                //!< end configuration (MDStep in atom::trajectory)
  atom **PermutationMap;      //!< gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
  DistanceMap **DistanceList; //!< distance list of each atom to each atom
  DistanceMap::iterator *StepList; //!< iterator to ascend through NearestNeighbours \a **DistanceList
  int *DoubleList;            //!< count of which sources want to move to this target, basically the injective measure (>1 -> not injective)
  DistanceMap::iterator *DistanceIterators; //!< marks which was the last picked target as injective candidate with smallest distance
  bool IsAngstroem;           //!< whether coordinates are in angstroem (true) or bohrradius (false)
  double *PenaltyConstants;   //!<  penalty constant in front of each term
};

#define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented 
enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output


/** The complete molecule.
 * Class incorporates number of types
 */
class molecule : public PointCloud , public Observable {
  friend molecule *NewMolecule();
  friend void DeleteMolecule(molecule *);
  public:
    const periodentafel * const elemente; //!< periodic table with each element
    atom *start;        //!< start of atom list
    atom *end;          //!< end of atom list
    bond *first;        //!< start of bond list
    bond *last;         //!< end of bond list
    int MDSteps;        //!< The number of MD steps in Trajectories
    int AtomCount;          //!< number of atoms, brought up-to-date by CountAtoms()
    int BondCount;          //!< number of atoms, brought up-to-date by CountBonds()
    int ElementCount;       //!< how many unique elements are therein
    int ElementsInMolecule[MAX_ELEMENTS]; //!< list whether element (sorted by atomic number) is alread present or not
    mutable int NoNonHydrogen;  //!< number of non-hydrogen atoms in molecule
    mutable int NoNonBonds;     //!< number of non-hydrogen bonds in molecule
    mutable int NoCyclicBonds;  //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
    double BondDistance;  //!< typical bond distance used in CreateAdjacencyList() and furtheron
    bool ActiveFlag;    //!< in a MoleculeListClass used to discern active from inactive molecules
    Vector Center;      //!< Center of molecule in a global box
    int IndexNr;        //!< index of molecule in a MoleculeListClass
    char name[MAXSTRINGSIZE];         //!< arbitrary name

  private:
    Cacheable<string> formula;
    moleculeId_t id;
  protected:
    molecule(const periodentafel * const teil);
    virtual ~molecule();


public:
  //getter and setter
  const std::string getName();
  moleculeId_t getId();
  void setId(moleculeId_t);
  void setName(const std::string);
  const std::string getFormula();
  std::string calcFormula();


  // re-definition of virtual functions from PointCloud
  const char * const GetName() const;
  Vector *GetCenter() const ;
  TesselPoint *GetPoint() const ;
  TesselPoint *GetTerminalPoint() const ;
  int GetMaxId() const;
  void GoToNext() const ;
  void GoToPrevious() const ;
  void GoToFirst() const ;
  void GoToLast() const ;
  bool IsEmpty() const ;
  bool IsEnd() const ;

  // templates for allowing global manipulation of all vectors
  template <typename res> void ActOnAllVectors( res (Vector::*f)() ) const;
  template <typename res> void ActOnAllVectors( res (Vector::*f)() const) const;
  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T), T t ) const;
  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T) const, T t ) const;
  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T&), T &t ) const;
  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T&) const, T &t ) const;
  template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const;
  template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const;
  template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v) const;
  template <typename res, typename T, typename U, typename V> void ActOnAllVectors( res (Vector::*f)(T, U, V) const, T t, U u, V v) const;

  // templates for allowing global manipulation of molecule with each atom as single argument
  template <typename res> void ActWithEachAtom( res (molecule::*f)(atom *) ) const;
  template <typename res> void ActWithEachAtom( res (molecule::*f)(atom *) const) const;

  // templates for allowing global copying of molecule with each atom as single argument
  template <typename res> void ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy) const;
  template <typename res> void ActOnCopyWithEachAtom( res (molecule::*f)(atom *) const, molecule *copy) const;

  // templates for allowing global manipulation of all atoms
  template <typename res, typename typ> void ActOnAllAtoms( res (typ::*f)() ) const;
  template <typename res, typename typ> void ActOnAllAtoms( res (typ::*f)() const) const;
  template <typename res, typename typ, typename T> void ActOnAllAtoms( res (typ::*f)(T), T t ) const;
  template <typename res, typename typ, typename T> void ActOnAllAtoms( res (typ::*f)(T) const, T t ) const;
  template <typename res, typename typ, typename T, typename U> void ActOnAllAtoms( res (typ::*f)(T, U), T t, U u ) const;
  template <typename res, typename typ, typename T, typename U> void ActOnAllAtoms( res (typ::*f)(T, U) const, T t, U u ) const;
  template <typename res, typename typ, typename T, typename U, typename V> void ActOnAllAtoms( res (typ::*f)(T, U, V), T t, U u, V v) const;
  template <typename res, typename typ, typename T, typename U, typename V> void ActOnAllAtoms( res (typ::*f)(T, U, V) const, T t, U u, V v) const;
  template <typename res, typename typ, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (typ::*f)(T, U, V, W), T t, U u, V v, W w) const;
  template <typename res, typename typ, typename T, typename U, typename V, typename W> void ActOnAllAtoms( res (typ::*f)(T, U, V, W) const, T t, U u, V v, W w) const;

  // templates for allowing conditional global copying of molecule with each atom as single argument
  template <typename res> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () ) const;
  template <typename res> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () const ) const;
  template <typename res> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) () ) const;
  template <typename res> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) () const ) const;
  template <typename res, typename T> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T), T t ) const;
  template <typename res, typename T> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T) const, T t ) const;
  template <typename res, typename T> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T), T t ) const;
  template <typename res, typename T> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T) const, T t ) const;
  template <typename res, typename T, typename U> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U), T t, U u ) const;
  template <typename res, typename T, typename U> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U) const, T t, U u ) const;
  template <typename res, typename T, typename U> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T, U), T t, U u ) const;
  template <typename res, typename T, typename U> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T, U) const, T t, U u ) const;
  template <typename res, typename T, typename U, typename V> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v ) const;
  template <typename res, typename T, typename U, typename V> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V) const, T t, U u, V v ) const;
  template <typename res, typename T, typename U, typename V> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v ) const;
  template <typename res, typename T, typename U, typename V> void ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) (T, U, V) const, T t, U u, V v ) const;

  // templates for allowing global manipulation of an array with one entry per atom
  void SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::* index) const;
  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::* index, void (*Setor)(T *, T *)) const;
  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::* index, void (*Setor)(T *, T *), T t) const;
  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::* index, void (*Setor)(T *, T *), T *t) const;
  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::* index, void (*Setor)(T *, T *)) const;
  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::* index, void (*Setor)(T *, T *), T t) const;
  template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int element::* index, void (*Setor)(T *, T *), T *t) const;
  template <typename T, typename typ> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ atom::*value) const;
  template <typename T, typename typ> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ atom::*value) const;
  template <typename T, typename typ> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ &vect ) const;
  template <typename T, typename typ> void SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ &vect ) const;

  // templates for allowing global manipulation of each atom by entries in an array
  template <typename T, typename typ, typename typ2> void SetAtomValueToIndexedArray ( T *array, int typ::*index, T typ2::*value ) const;
  template <typename T, typename typ> void SetAtomValueToValue ( T value, T typ::*ptr ) const;

  template <typename res, typename typ> res SumPerAtom(res (typ::*f)() ) const;
  template <typename res, typename typ> res SumPerAtom(res (typ::*f)() const ) const;
  template <typename res, typename typ, typename T> res SumPerAtom(res (typ::*f)(T) , T t ) const;
  template <typename res, typename typ, typename T> res SumPerAtom(res (typ::*f)(T) const, T t ) const;

  /// remove atoms from molecule.
  bool AddAtom(atom *pointer);
  bool RemoveAtom(atom *pointer);
  bool UnlinkAtom(atom *pointer);
  bool CleanupMolecule();

  /// Add/remove atoms to/from molecule.
  atom * AddCopyAtom(atom *pointer);
  bool AddXYZFile(string filename);
  bool AddHydrogenReplacementAtom(bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem);
  bond * AddBond(atom *first, atom *second, int degree = 1);
  bool RemoveBond(bond *pointer);
  bool RemoveBonds(atom *BondPartner);

  /// Find atoms.
  atom * FindAtom(int Nr) const;
  atom * AskAtom(string text);

  /// Count and change present atoms' coordination.
  void CountAtoms();
  void CountElements();
  void CalculateOrbitals(class config &configuration);
  bool CenterInBox();
  bool BoundInBox();
  void CenterEdge(Vector *max);
  void CenterOrigin();
  void CenterPeriodic();
  void CenterAtVector(Vector *newcenter);
  void Translate(const Vector *x);
  void TranslatePeriodically(const Vector *trans);
  void Mirror(const Vector *x);
  void Align(Vector *n);
  void Scale(const double ** const factor);
  void DeterminePeriodicCenter(Vector &center);
  Vector * DetermineCenterOfGravity();
  Vector * DetermineCenterOfAll() const;
  void SetNameFromFilename(const char *filename);
  void SetBoxDimension(Vector *dim);
  void ScanForPeriodicCorrection();
  bool VerletForceIntegration(char *file, config &configuration);
  void Thermostats(config &configuration, double ActualTemp, int Thermostat);
  void PrincipalAxisSystem(bool DoRotate);
  double VolumeOfConvexEnvelope(bool IsAngstroem);

  double ConstrainedPotential(struct EvaluatePotential &Params);
  double MinimiseConstrainedPotential(atom **&permutation, int startstep, int endstep, bool IsAngstroem);
  void EvaluateConstrainedForces(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
  bool LinearInterpolationBetweenConfiguration(int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity);
        
  bool CheckBounds(const Vector *x) const;
  void GetAlignvector(struct lsq_params * par) const;

  /// Initialising routines in fragmentation
  void CreateAdjacencyListFromDbondFile(ifstream *output);
  void CreateAdjacencyList(double bonddistance, bool IsAngstroem, void (BondGraph::*f)(BondedParticle * const , BondedParticle * const , double &, double &, bool), BondGraph *BG = NULL);
  int CorrectBondDegree() const;
  void OutputBondsList() const;
  void CyclicBondAnalysis() const;
  void OutputGraphInfoPerAtom() const;
  void OutputGraphInfoPerBond() const;


  // Graph analysis
  MoleculeLeafClass * DepthFirstSearchAnalysis(class StackClass<bond *> *&BackEdgeStack) const;
  void CyclicStructureAnalysis(class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize) const;
  bool PickLocalBackEdges(atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack) const;
  bond * FindNextUnused(atom *vertex) const;
  void SetNextComponentNumber(atom *vertex, int nr) const;
  void ResetAllBondsToUnused() const;
  int CountCyclicBonds();
  bool CheckForConnectedSubgraph(KeySet *Fragment);
  string GetColor(enum Shading color) const;
  bond * CopyBond(atom *left, atom *right, bond *CopyBond);


  molecule *CopyMolecule();
  molecule* CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const;

  /// Fragment molecule by two different approaches:
  int FragmentMolecule(int Order, config *configuration);
  bool CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
  bool StoreBondsToFile(char *path, char *filename);
  bool StoreAdjacencyToFile(char *path, char *filename);
  bool CheckAdjacencyFileAgainstMolecule(char *path, atom **ListOfAtoms);
  bool ParseOrderAtSiteFromFile(char *path);
  bool StoreOrderAtSiteFile(char *path);
  bool StoreForcesFile(MoleculeListClass *BondFragments, char *path, int *SortIndex);
  bool CreateMappingLabelsToConfigSequence(int *&SortIndex);
  void BreadthFirstSearchAdd(molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
  /// -# BOSSANOVA
  void FragmentBOSSANOVA(Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize);
  int PowerSetGenerator(int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet);
  bool BuildInducedSubgraph(const molecule *Father);
  molecule * StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem);
  void SPFragmentGenerator(struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder);
  int LookForRemovalCandidate(KeySet *&Leaf, int *&ShortestPathList);
  int GuesstimateFragmentCount(int order);

  // Recognize doubly appearing molecules in a list of them
  int * IsEqualToWithinThreshold(molecule *OtherMolecule, double threshold);
  int * GetFatherSonAtomicMap(molecule *OtherMolecule);

  // Output routines.
  bool Output(ofstream * const output);
  bool OutputTrajectories(ofstream * const output);
  void OutputListOfBonds() const;
  bool OutputXYZ(ofstream * const output) const;
  bool OutputTrajectoriesXYZ(ofstream * const output);
  bool Checkout(ofstream * const output) const;
  bool OutputTemperatureFromTrajectories(ofstream * const output, int startstep, int endstep);

  // Manipulation routines
  void flipActiveFlag();

  private:
  int last_atom;      //!< number given to last atom
  mutable atom *InternalPointer;  //!< internal pointer for PointCloud
};

molecule *NewMolecule();
void DeleteMolecule(molecule* mol);

#include "molecule_template.hpp"

/** A list of \a molecule classes.
 */
class MoleculeListClass : public Observable {
  public:
    MoleculeList ListOfMolecules; //!< List of the contained molecules
    int MaxIndex;

  MoleculeListClass(World *world);
  ~MoleculeListClass();

  bool AddHydrogenCorrection(char *path);
  bool StoreForcesFile(char *path, int *SortIndex);
  void insert(molecule *mol);
  molecule * ReturnIndex(int index);
  bool OutputConfigForListOfFragments(config *configuration, int *SortIndex);
  int NumberOfActiveMolecules();
  void Enumerate(ostream *out);
  void Output(ofstream *out);
  void DissectMoleculeIntoConnectedSubgraphs(const periodentafel * const periode, config * const configuration);
  int CountAllAtoms() const;

  // Methods moved here from the menus
  // TODO: more refactoring needed on these methods
  void flipChosen();
  void createNewMolecule(periodentafel *periode);
  void loadFromXYZ(periodentafel *periode);
  void setMoleculeFilename();
  void parseXYZIntoMolecule();
  void eraseMolecule();


  // merging of molecules
  bool SimpleMerge(molecule *mol, molecule *srcmol);
  bool SimpleAdd(molecule *mol, molecule *srcmol);
  bool SimpleMultiMerge(molecule *mol, int *src, int N);
  bool SimpleMultiAdd(molecule *mol, int *src, int N);
  bool ScatterMerge(molecule *mol, int *src, int N);
  bool EmbedMerge(molecule *mol, molecule *srcmol);

  private:
  World *world; //!< The world this List belongs to. Needed to avoid deadlocks in the destructor
};


/** A leaf for a tree of \a molecule class
 * Wraps molecules in a tree structure
 */
class MoleculeLeafClass {
  public:
    molecule *Leaf;                   //!< molecule of this leaf
    //MoleculeLeafClass *UpLeaf;        //!< Leaf one level up
    //MoleculeLeafClass *DownLeaf;      //!< First leaf one level down
    MoleculeLeafClass *previous;  //!< Previous leaf on this level
    MoleculeLeafClass *next;      //!< Next leaf on this level

  //MoleculeLeafClass(MoleculeLeafClass *Up, MoleculeLeafClass *Previous);
  MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf);
  ~MoleculeLeafClass();

  bool AddLeaf(molecule *ptr, MoleculeLeafClass *Previous);
  bool FillBondStructureFromReference(const molecule * const reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList = false);
  bool FillRootStackForSubgraphs(KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter);
  bool AssignKeySetsToFragment(molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList = false);
  bool FillListOfLocalAtoms(atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList);
  void TranslateIndicesToGlobalIDs(Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph);
  int Count() const;
};


#endif /*MOLECULES_HPP_*/