source: src/molecule.hpp@ 47f0e4

/** \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 <boost/bimap/bimap.hpp>
#include <boost/bimap/unordered_set_of.hpp>
#include <boost/bimap/multiset_of.hpp>
#include <boost/optional.hpp>
#include <boost/shared_ptr.hpp>

#include "AtomIdSet.hpp"
#include "Atom/AtomSet.hpp"
#include "CodePatterns/Cacheable.hpp"
#include "CodePatterns/Observer/Observable.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "Fragmentation/HydrogenSaturation_enum.hpp"
#include "Formula.hpp"
#include "Helpers/defs.hpp"
#include "IdPool_policy.hpp"
#include "IdPool.hpp"
#include "Shapes/Shape.hpp"
#include "types.hpp"

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

class atom;
class bond;
class BondedParticle;
class BondGraph;
class DepthFirstSearchAnalysis;
class element;
class ForceMatrix;
class Graph;
class LinkedCell_deprecated;
class ListOfLocalAtoms_t;
class molecule;
class MoleculeLeafClass;
class MoleculeUnittest;
class RealSpaceMatrix;
class Vector;

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

/** External function to remove all atoms since this will also delete the molecule
 *
 * \param _mol ref pointer to molecule to destroy
 */
void removeAtomsinMolecule(molecule *&_mol);

/** The complete molecule.
 * Class incorporates number of types
 */
class molecule : public Observable
{
  //!> grant unit test access
  friend class MoleculeUnittest;
  //!> function may access cstor
  friend molecule *NewMolecule();
  //!> function may access dstor
  friend void DeleteMolecule(molecule *);

public:
  typedef AtomIdSet::atomIdSet atomIdSet;
  typedef AtomIdSet::iterator iterator;
  typedef AtomIdSet::const_iterator const_iterator;

  int MDSteps; //!< The number of MD steps in Trajectories
  mutable int NoNonBonds; //!< number of non-hydrogen bonds in molecule
  mutable int NoCyclicBonds; //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
  bool ActiveFlag; //!< in a MoleculeListClass used to discern active from inactive molecules
  int IndexNr; //!< index of molecule in a MoleculeListClass
  char name[MAXSTRINGSIZE]; //!< arbitrary name

private:
  Formula formula;
  size_t NoNonHydrogen; //!< number of non-hydrogen atoms in molecule
  int BondCount; //!< number of atoms, brought up-to-date by doCountBonds()
  moleculeId_t id;
  AtomIdSet atomIds; //<!set of atomic ids to check uniqueness of atoms
  IdPool<atomId_t, uniqueId> atomIdPool;  //!< pool of internal ids such that way may guarantee uniqueness
  typedef std::map<atomId_t,atom *> LocalToGlobalId_t;
  LocalToGlobalId_t LocalToGlobalId; //!< internal map to ease FindAtom

protected:

  molecule();
  virtual ~molecule();

public:

  /******* Notifications *******/

  //!> enumeration of present notification types: only insertion/removal of atoms or molecules
  enum NotificationType {
    AtomInserted,
    AtomRemoved,
    AtomNrChanged,
    AtomMoved,
    FormulaChanged,
    MoleculeCenterChanged,
    MoleculeNameChanged,
    IndexChanged,
    BoundingBoxChanged,
    AboutToBeRemoved,
    SelectionChanged,
    NotificationType_MAX
  };

  //>! access to last changed element (atom)
  const atomId_t lastChangedAtomId() const
  { return _lastchangedatomid; }

public:
  //getter and setter
  const std::string getName() const;
  int getAtomCount() const;
  size_t getNoNonHydrogen() const{
    return NoNonHydrogen;
  }

  int getBondCount() const{
    return BondCount;
  }
  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);

  World::AtomComposite getAtomSet();
  World::ConstAtomComposite getAtomSet() const;

  // simply pass on all functions to AtomIdSet
  iterator begin() {
    return atomIds.begin();
  }
  const_iterator begin() const
  {
    return atomIds.begin();
  }
  iterator end()
  {
    return atomIds.end();
  }
  const_iterator end() const
  {
    return atomIds.end();
  }
  bool empty() const
  {
    return atomIds.empty();
  }
  size_t size() const
  {
    return atomIds.size();
  }
  const_iterator find(atom * key) const
  {
    return atomIds.find(key);
  }

  /** Returns the set of atomic ids contained in this molecule.
   *
   * @return set of atomic ids
   */
  const atomIdSet & getAtomIds() const {
    return atomIds.getAtomIds();
  }

  std::pair<iterator, bool> insert(atom * const key);

  /** Predicate whether given \a key is contained in this molecule.
   *
   * @param key atom to check
   * @return true - is contained, false - else
   */
  bool containsAtom(const atom* key) const
  {
    return atomIds.contains(key);
  }

  /** Predicate whether given \a id is contained in this molecule.
   *
   * @param id atomic id to check
   * @return true - is contained, false - else
   */
  bool containsAtom(const atomId_t id) const
  {
    return atomIds.contains(id);
  }

private:
  friend void atom::removeFromMolecule();
  /** Erase an atom from the list.
   * \note This should only be called by atom::removeFromMolecule(),
   * otherwise it is not assured that the atom knows about it.
   *
   * @param loc locator to atom in list
   * @return iterator to just after removed item (compliant with standard)
   */
  const_iterator erase(const_iterator loc);

  /** Erase an atom from the list.
   * \note This should only be called by atom::removeFromMolecule(),
   * otherwise it is not assured that the atom knows about it.
   *
   * @param *key key to atom in list
   * @return iterator to just after removed item (compliant with standard)
   */
  const_iterator erase(atom * key);

private:
  friend bool atom::changeNr(int newId);
  /**
   * used when changing an ParticleInfo::Nr.
   * Note that this number is local with this molecule.
   * Unless you are calling this method from inside an atom don't fiddle with the third parameter.
   *
   * @param oldNr old Nr
   * @param newNr new Nr to set
   * @param *target ref to atom
   * @return indicates wether the change could be done or not.
   */
  bool changeAtomNr(int oldNr, int newNr, atom* target=0);

  friend bool atom::changeId(atomId_t newId);
  /**
   * used when changing an ParticleInfo::Id.
   * Note that this number is global (and the molecule uses it to know which atoms belong to it)
   *
   * @param oldId old Id
   * @param newId new Id to set
   * @return indicates wether the change could be done or not.
   */
  bool changeAtomId(int oldId, int newId);

  /** Updates the internal lookup fro local to global indices.
   *
   * \param pointer pointer to atom
   */
  void InsertLocalToGlobalId(atom * const pointer);

  /** Sets the name of the atom.
   *
   * The name is set via its element symbol and its internal ParticleInfo::Nr.
   *
   * @param _atom atom whose name to set
   */
  void setAtomName(atom *_atom) const;

  //!> grant World (only) access to selection state changers
  friend class World;

  /** Sets the internal selection state to true.
   *
   */
  void select();

  /** Unsets the internal selection state to true.
   *
   */
  void unselect();

public:

  /** Getter to internal selection status.
   *
   * \return true - molecule is selected, false - else
   */
  bool getSelected() const { return selected; }

  /** Structure for the required information on the bounding box.
   *
   */
  struct BoundingBoxInfo {
    //!> position of center
    Vector position;
    //!> radius of sphere
    double radius;

    /** Equivalence operator for bounding box.
     *
     * \return true - both bounding boxes have same position and radius
     */
    bool operator==(const BoundingBoxInfo &_other) const
    {  return (radius == _other.radius) && (position == _other.position); }

    /** Inequivalence operator for bounding box.
     *
     * \return true - bounding boxes have either different positions or different radii or both
     */
    bool operator!=(const BoundingBoxInfo &_other) const
    { return !(*this == _other); }
  };

private:

  /** Returns the current bounding box.
   *
   * \return Shape with center and extension of box
   */
  BoundingBoxInfo updateBoundingBox() const;

  /** Returns the current center of the molecule.
   *
   * \return center
   */
  Vector updateMoleculeCenter() const;

  // stuff for keeping bounding box up-to-date efficiently

  //!> Cacheable for the bounding box, ptr such that
  boost::shared_ptr< Cacheable<BoundingBoxInfo> > BoundingBox;

  //!> Cacheable for the bounding box, ptr such that
  boost::shared_ptr< Cacheable<Vector> > MoleculeCenter;

  /** Bimap storing atomic ids and the component per axis.
   *
   * We need a bimap in order to have the components sorted and be able to
   * access max and min values in linear time and also access the ids in
   * constant time in order to update the map, when atoms move, are inserted,
   * or removed.
   */
  typedef boost::bimaps::bimap<
          boost::bimaps::set_of< atomId_t >,
          boost::bimaps::multiset_of< double, std::greater<double> >
      > AtomDistanceMap_t;
  std::vector<AtomDistanceMap_t> BoundingBoxSweepingAxis;

  //!> typedef for a map with current bond counts per atom
  typedef std::map<atomId_t, size_t> BondCountsPerAtom_t;

  //!> current bond counts per atom to update the BondCount
  BondCountsPerAtom_t BondCountsPerAtom;

  //!> typedef for a map with current element per atom
  typedef std::map<atomId_t, atomicNumber_t> ElementPerAtom_t;

  //!> current element per atom to update the BondCount
  ElementPerAtom_t ElementPerAtom;

  //!> make setMolecule friend to access associateAtomWithMolecule()
  friend void atom::setMolecule(molecule *);

  /** Helper function only to be called by specific atom function.
   *
   * \param _atom atom to be added to this molecule
   */
  void associateAtomWithMolecule(atom *_atom);

  /** Helper function only to be called by specific atom function.
   *
   * \param _atom atom to be added to this molecule
   */
  void disassociateAtomWithMolecule(atom *_atom);

public:

  /** Returns the current bounding box of this molecule.
   *
   * \return bounding box info with center and radius
   */
  BoundingBoxInfo getBoundingBox() const;

  /** Function to create a bounding spherical shape for the currently associated atoms.
   *
   * \param boundary extra boundary of shape around (i.e. distance between outermost atom
   *        and the shape's surface)
   */
  Shape getBoundingSphere(const double boundary = 0.) const;

  /** Creates the bounding box by adding van der Waals-Spheres around every atom.
   *
   * \param scale extra scale parameter to enlarge the spheres artifically
   */
  Shape getBoundingShape(const double scale = 1.) const;

  /** Returns the current center of this molecule.
   *
   * \return center of the molecule
   */
  Vector getMoleculeCenter() 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 AddHydrogenReplacementAtom(bond::ptr Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem);
  bond::ptr AddBond(atom *first, atom *second, int degree = 1);
  bool hasBondStructure() const;

  /// Find atoms.
  atom * FindAtom(int Nr) const;
  atom * AskAtom(std::string text);
  bool isInMolecule(const atom * const _atom) const;

  /// Count and change present atoms' coordination.
  bool CenterInBox();
  bool BoundInBox();
  void CenterEdge();
  void CenterOrigin();
  void CenterPeriodic();
  void CenterAtVector(const Vector &newcenter);
  void Translate(const Vector &x);
  void TranslatePeriodically(const Vector &trans);
  void Mirror(const Vector &x);
  void Align(const Vector &n);
  void Scale(const double *factor);
  void DeterminePeriodicCenter(Vector &center, const enum HydrogenTreatment _treatment = ExcludeHydrogen);
  const Vector DetermineCenterOfGravity() const;
  const Vector DetermineCenterOfAll() const;
  void SetNameFromFilename(const char *filename);
  bool ScanForPeriodicCorrection();
  double VolumeOfConvexEnvelope(bool IsAngstroem);
  RealSpaceMatrix getInertiaTensor() const;
  void RotateToPrincipalAxisSystem(const Vector &Axis);

  bool CheckBounds(const Vector *x) const;
  void GetAlignvector(struct lsq_params * par) const;

  /// Initialising routines in fragmentation
  void OutputBondsList() const;

  bond::ptr CopyBond(atom *left, atom *right, bond::ptr CopyBond);

  molecule *CopyMolecule(const Vector &offset = zeroVec);
  molecule* CopyMoleculeFromSubRegion(const Shape&);

  /// Fragment molecule by two different approaches:
  bool StoreBondsToFile(std::string filename, std::string path = "");
  bool CreateFatherLookupTable(ListOfLocalAtoms_t &LookupTable, int count = 0);

  // Recognize doubly appearing molecules in a list of them
  int * GetFatherSonAtomicMap(const molecule * const OtherMolecule);
  bool FillBondStructureFromReference(const molecule * const reference, ListOfLocalAtoms_t &ListOfLocalAtoms, bool FreeList = false);
  bool FillListOfLocalAtoms(ListOfLocalAtoms_t &ListOfLocalAtoms, const int GlobalAtomCount);

  // Output routines.
  bool Output(std::ostream * const output) const;
  void OutputListOfBonds() const;

  // Manipulation routines
  void flipActiveFlag();

  virtual void update(Observable *publisher);
  virtual void recieveNotification(Observable *publisher, Notification_ptr notification);
  virtual void subjectKilled(Observable *publisher);

private:
  //!> id of last atom that signalled changed associated with this molecule
  atomId_t _lastchangedatomid;

  int last_atom; //!< number given to last atom

  //!> center of the molecule
  Vector molcenter;

  //!> internal state whether atom is selected or not
  bool selected;
};

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



#endif /*MOLECULES_HPP_*/