source: src/molecule.hpp@ a0064e

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

#ifndef MOLECULES_HPP_
#define MOLECULES_HPP_

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

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

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

#include <string>

#include "types.hpp"
#include "graph.hpp"
#include "PointCloud.hpp"
#include "Patterns/Observer.hpp"
#include "Patterns/ObservedIterator.hpp"
#include "Patterns/Cacheable.hpp"
#include "Formula.hpp"
#include "AtomSet.hpp"

#include "Descriptors/MoleculeDescriptor_impl.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;
class Shape;

/******************************** 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
};

/** The complete molecule.
 * Class incorporates number of types
 */
class molecule : public PointCloud, public Observable
{
  friend molecule *NewMolecule();
  friend void DeleteMolecule(molecule *);

public:
  typedef ATOMSET(std::list) atomSet;
  typedef std::set<atomId_t> atomIdSet;
  typedef ObservedIterator<atomSet> iterator;
  typedef atomSet::const_iterator const_iterator;

  const periodentafel * const elemente; //!< periodic table with each element
  // old deprecated atom handling
  //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()
  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:
  Formula formula;
  Cacheable<int> AtomCount;
  moleculeId_t id;
  atomSet atoms; //<!list of atoms
  atomIdSet atomIds; //<!set of atomic ids to check uniqueness of atoms
protected:
  //void CountAtoms();
  /**
   * this iterator type should be used for internal variables, \
     * since it will not lock
   */
  typedef atomSet::iterator internal_iterator;

  molecule(const periodentafel * const teil);
  virtual ~molecule();

public:
  //getter and setter
  const std::string getName() const;
  int getAtomCount() const;
  int doCountAtoms();
  moleculeId_t getId() const;
  void setId(moleculeId_t);
  void setName(const std::string);
  const Formula &getFormula() const;
  unsigned int getElementCount() const;
  bool hasElement(const element*) const;
  bool hasElement(atomicNumber_t) const;
  bool hasElement(const std::string&) const;

  virtual bool changeId(atomId_t newId);

  TesselPoint * getValue(const_iterator &rhs) const;
  TesselPoint * getValue(iterator &rhs) const;
  iterator begin();
  const_iterator begin() const;
  iterator end();
  const_iterator end() const;
  bool empty() const;
  size_t size() const;
  const_iterator erase(const_iterator loc);
  const_iterator erase(atom * key);
  const_iterator find(atom * key) const;
  pair<iterator, bool> insert(atom * const key);
  bool containsAtom(atom* key);

  // re-definition of virtual functions from PointCloud
  const char * const GetName() const;
  Vector *GetCenter() const;
  TesselPoint *GetPoint() const;
  int GetMaxId() const;
  void GoToNext() const;
  void GoToFirst() const;
  bool IsEmpty() const;
  bool IsEnd() 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);
  bool hasBondStructure() const;
  unsigned int CountBonds() const;

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

  /// Count and change present atoms' coordination.
  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() const;
  Vector * DetermineCenterOfAll() const;
  Vector * DetermineCenterOfBox() const;
  void SetNameFromFilename(const char *filename);
  void SetBoxDimension(Vector *dim);
  void ScanForPeriodicCorrection();
  bool VerletForceIntegration(char *file, config &configuration, const size_t offset);
  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, std::string 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(std::deque<bond *> *&BackEdgeStack) const;
  void CyclicStructureAnalysis(std::deque<bond *> *BackEdgeStack, int *&MinimumRingSize) const;
  bool PickLocalBackEdges(atom **ListOfLocalAtoms, std::deque<bond *> *&ReferenceStack, std::deque<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() const;
  molecule* CopyMoleculeFromSubRegion(const Shape&) const;

  /// Fragment molecule by two different approaches:
  int FragmentMolecule(int Order, std::string &prefix);
  bool CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, std::string path = "");
  bool StoreBondsToFile(std::string filename, std::string path = "");
  bool StoreAdjacencyToFile(std::string filename, std::string path = "");
  bool CheckAdjacencyFileAgainstMolecule(std::string &path, atom **ListOfAtoms);
  bool ParseOrderAtSiteFromFile(std::string &path);
  bool StoreOrderAtSiteFile(std::string &path);
  bool StoreForcesFile(MoleculeListClass *BondFragments, std::string &path, int *SortIndex);
  bool CreateMappingLabelsToConfigSequence(int *&SortIndex);
  bool CreateFatherLookupTable(atom **&LookupTable, int count = 0);
  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 * GetFatherSonAtomicMap(molecule *OtherMolecule);

  // Output routines.
  bool Output(std::ostream * const output) const;
  bool OutputTrajectories(ofstream * const output) const;
  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:
  void init_DFS(struct DFSAccounting&) const;
  int last_atom; //!< number given to last atom
  mutable internal_iterator InternalPointer; //!< internal pointer for PointCloud
};

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

/** 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(std::string &path);
  bool StoreForcesFile(std::string &path, int *SortIndex);
  void insert(molecule *mol);
  void erase(molecule *mol);
  molecule * ReturnIndex(int index);
  bool OutputConfigForListOfFragments(std::string &prefix, int *SortIndex);
  int NumberOfActiveMolecules();
  void Enumerate(ostream *out);
  void Output(ofstream *out);
  int CountAllAtoms() const;

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

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, 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 GlobalAtomCount, bool &FreeList);
  void TranslateIndicesToGlobalIDs(Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph);
  int Count() const;
};

#endif /*MOLECULES_HPP_*/