Diff [2a03b00b33c43f417fdf3452727473a2043f2638:ee4f2da6990041aedcb8b339896459f5b4e9b05b] for / – MoleCuilder

LinearAlgebra/src/LinearAlgebra/Makefile.am

-              r2a03b0
+              ree4f2d
 INCLUDES = -I$(top_srcdir)/src
 AM_LDFLAGS = ${CodePatterns_LIBS} -ldl
+AM_LDFLAGS = -ldl
 AM_CPPFLAGS = $(BOOST_CPPFLAGS) ${CodePatterns_CFLAGS}
 …
 libLinearAlgebra_la_LIBADD = \
         $(BOOST_SERIALIZATION_LDFLAGS) $(BOOST_SERIALIZATION_LIBS) \
+        ${CodePatterns_LIBS} \
         $(GSL_LIBS)
 nobase_libLinearAlgebra_la_include_HEADERS = ${LINALGHEADER}

LinearAlgebra/src/LinearAlgebra/Plane.cpp

-              r2a03b0
+              ree4f2d
   Vector x1 = y1 - y2;
   Vector x2 = y3 - y2;
   if ((x1.Norm() <= LINALG_MYEPSILON())) {
+  if ((x1.NormSquared() <= LINALG_MYEPSILON())) {
     throw LinearDependenceException() << LinearAlgebraVectorPair( make_pair(&y1, &y2) );
+  }
   if ((x2.Norm() <= LINALG_MYEPSILON())) {
+  if ((x2.NormSquared() <= LINALG_MYEPSILON())) {
     throw LinearDependenceException() << LinearAlgebraVectorPair( make_pair(&y2, &y3) );
+  }
+  if((fabs(x1.Angle(x2)) <= LINALG_MYEPSILON())) {
+  const Vector lineardependent = x1.Projection(x2) - x1;
+  if((lineardependent.NormSquared() <= LINALG_MYEPSILON())) {
     throw LinearDependenceException() << LinearAlgebraVectorPair( make_pair(&x1, &x2) );
+  }

LinearAlgebra/src/LinearAlgebra/Vector.cpp

-              r2a03b0
+              ree4f2d
 /** Checks whether vector has all components zero.
+ * @param epsilon numerical tolerance for equality
  * @return true - vector is zero, false - vector is not
  */
 bool Vector::IsZero() const
+{
   return (fabs(at(0))+fabs(at(1))+fabs(at(2)) < LINALG_MYEPSILON());
+bool Vector::IsZero(const double _epsilon) const
+{
+  return (fabs(at(0))+fabs(at(1))+fabs(at(2)) < _epsilon);
 };
 /** Checks whether vector has length of 1.
+ * @param epsilon numerical tolerance for equality
  * @return true - vector is normalized, false - vector is not
  */
 bool Vector::IsOne() const
+{
   return (fabs(Norm() - 1.) < LINALG_MYEPSILON());
+bool Vector::IsOne(const double _epsilon) const
+{
+  return (fabs(Norm() - 1.) < _epsilon);
 };
 /** Checks whether vector is normal to \a *normal.
+ * @param normal other supposedly normal vector
+ * @param epsilon numerical tolerance for equality
  * @return true - vector is normalized, false - vector is not
  */
 bool Vector::IsNormalTo(const Vector &normal) const
+{
   if (ScalarProduct(normal) < LINALG_MYEPSILON())
+bool Vector::IsNormalTo(const Vector &normal, const double _epsilon) const
+{
+  if (ScalarProduct(normal) < _epsilon)
     return true;
   else
 …
 };
+/** Checks whether vector is parallel to \a *parallel.
+ * @param parallel other supposedly parallel vector
+ * @param epsilon numerical tolerance for equality
+ * @return true - vector is parallel/linear dependent, false - vector is not
+ */
+bool Vector::IsParallelTo(const Vector &parallel, const double _epsilon) const
+{
+  if (1.-fabs(ScalarProduct(parallel)/parallel.Norm()/this->Norm()) < _epsilon)
+    return true;
+  else
+    return false;
+};
 /** Checks whether vector is normal to \a *normal.
+ * @param a other vector
+ * @param epsilon numerical tolerance for equality
  * @return true - vector is normalized, false - vector is not
  */
 bool Vector::IsEqualTo(const Vector &a) const
+bool Vector::IsEqualTo(const Vector &a, const double _epsilon) const
+{
   bool status = true;
   for (int i=0;i<NDIM;i++) {
     if (fabs(at(i) - a[i]) > LINALG_MYEPSILON())
+    if (fabs(at(i) - a[i]) > _epsilon)
       status = false;
+  }

LinearAlgebra/src/LinearAlgebra/Vector.hpp

-              r2a03b0
+              ree4f2d
   double ScalarProduct(const Vector &y) const;
   double Angle(const Vector &y) const;
+  bool IsZero() const;
+  bool IsOne() const;
+  bool IsNormalTo(const Vector &normal) const;
+  bool IsEqualTo(const Vector &a) const;
+  bool IsZero(const double epsilon = LINALG_MYEPSILON()) const;
+  bool IsOne(const double epsilon = LINALG_MYEPSILON()) const;
+  bool IsNormalTo(const Vector &normal, const double epsilon = LINALG_MYEPSILON()) const;
+  bool IsParallelTo(const Vector &normal, const double epsilon = LINALG_MYEPSILON()) const;
+  bool IsEqualTo(const Vector &a, const double epsilon = LINALG_MYEPSILON()) const;
   void AddVector(const Vector &y);

LinearAlgebra/src/LinearAlgebra/VectorContent.cpp

-              r2a03b0
+              ree4f2d
+{
   double factor = Norm();
+  (*this) *= 1/factor;
+  if (fabs(factor) > LINALG_MYEPSILON())
+    (*this) *= 1./factor;
 };

Makefile.am

-              r2a03b0
+              ree4f2d
 doc:
         mkdir -p ${DX_DOCDIR}
         cd doc && make doxygen-doc
+        cd doc && $(MAKE) doxygen-doc
 extracheck:
         cd tests && make extracheck
+        cd tests && $(MAKE) extracheck
 installextracheck:
         cd tests && make installextracheck
+        cd tests && $(MAKE) installextracheck
 unity:
         cd src && make unity
+        cd src && $(MAKE) unity

src/Actions/FragmentationAction/FragmentationAction.cpp

-              r2a03b0
+              ree4f2d
 #include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp"
 #include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"
+#include "Fragmentation/Exporters/SaturatedBond.hpp"
+#include "Fragmentation/Exporters/SaturatedFragment.hpp"
+#include "Fragmentation/Exporters/SaturationDistanceMaximizer.hpp"
 #include "Fragmentation/Fragmentation.hpp"
 #include "Fragmentation/Graph.hpp"
 …
+  }
+  // create global saturation positions map
+  SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
+  {
+    // go through each atom
+    for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
+        iter != world.endAtomSelection(); ++iter) {
+      const atom * const _atom = iter->second;
+      // skip hydrogens if treated special
+      const enum HydrogenTreatment treatment =  params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
+      if ((treatment == ExcludeHydrogen) && (_atom->getType()->getAtomicNumber() == 1)) {
+        LOG(4, "DEBUG: Skipping hydrogen atom " << *_atom);
+        continue;
+      }
+      // get the valence
+      unsigned int NumberOfPoints = _atom->getElement().getNoValenceOrbitals();
+      LOG(3, "DEBUG: There are " << NumberOfPoints
+          << " places to fill in in total for this atom " << *_atom << ".");
+      // check whether there are any bonds with degree larger than 1
+      unsigned int SumOfDegrees = 0;
+      bool PresentHigherBonds = false;
+      const BondList &bondlist = _atom->getListOfBonds();
+      for (BondList::const_iterator bonditer = bondlist.begin();
+          bonditer != bondlist.end(); ++bonditer) {
+        SumOfDegrees += (*bonditer)->getDegree();
+        PresentHigherBonds |= (*bonditer)->getDegree() > 1;
+      }
+      // check whether there are alphas to maximize the hydrogens distances
+      SaturationDistanceMaximizer::position_bins_t position_bins;
+      {
+        // gather all bonds and convert to SaturatedBonds
+        SaturationDistanceMaximizer::PositionContainers_t CutBonds;
+        for (BondList::const_iterator bonditer = bondlist.begin();
+            bonditer != bondlist.end(); ++bonditer) {
+          CutBonds.push_back(
+              SaturatedBond::ptr(new SaturatedBond(*(bonditer->get()), *_atom) )
+            );
+        }
+        SaturationDistanceMaximizer maximizer(CutBonds);
+        if (PresentHigherBonds) {
+          // then find best alphas
+          maximizer();
+        } else {
+          // if no higher order bonds, we simply gather the scaled positions
+        }
+        position_bins = maximizer.getAllPositionBins();
+        LOG(4, "DEBUG: Positions for atom " << *_atom << " are " << position_bins);
+      }
+      // convert into the desired entry in the map
+      SaturatedFragment::SaturationsPositionsPerNeighbor_t positions_per_neighbor;
+      {
+        BondList::const_iterator bonditer = bondlist.begin();
+        SaturationDistanceMaximizer::position_bins_t::const_iterator biniter =
+            position_bins.begin();
+        for (;bonditer != bondlist.end(); ++bonditer, ++biniter) {
+          const atom * const OtherAtom = (*bonditer)->GetOtherAtom(_atom);
+          std::pair<
+              SaturatedFragment::SaturationsPositionsPerNeighbor_t::iterator,
+              bool
+              > inserter;
+          // check whether we treat hydrogen special
+          if ((treatment == ExcludeHydrogen) && (OtherAtom->getType()->getAtomicNumber() == 1)) {
+            // if hydrogen, forget rescaled position and use original one
+            inserter =
+                positions_per_neighbor.insert(
+                    std::make_pair(
+                        OtherAtom->getId(),
+                        SaturatedFragment::SaturationsPositions_t(
+, OtherAtom->getPosition() - _atom->getPosition())
+                    )
+                );
+          } else {
+            inserter =
+                positions_per_neighbor.insert(
+                    std::make_pair(
+                        OtherAtom->getId(),
+                        SaturatedFragment::SaturationsPositions_t(
+                            biniter->begin(),
+                            biniter->end())
+                    )
+                );
+          }
+          // if already pressent, add to this present list
+          ASSERT (inserter.second,
+              "FragmentationAction::performCall() - other atom "
+              +toString(*OtherAtom)+" already present?");
+        }
+        // bonditer follows nicely
+        ASSERT( biniter == position_bins.end(),
+            "FragmentationAction::performCall() - biniter is out of step, it still points at bond "
+            +toString(*biniter)+".");
+      }
+      // and insert
+      globalsaturationpositions.insert(
+          std::make_pair( _atom->getId(),
+              positions_per_neighbor
+          ));
+    }
+  }
+  {
     const enum HydrogenSaturation saturation =  params.DoSaturation.get() ? DoSaturate : DontSaturate;
 …
     if (params.types.get().size() != 0) {
       // store molecule's fragment to file
       ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation);
+      ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
       exporter.setPrefix(params.prefix.get());
       exporter.setOutputTypes(params.types.get());
 …
     } else {
       // store molecule's fragment in FragmentJobQueue
       ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation);
+      ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
       exporter.setLevel(params.level.get());
       exporter();

src/Actions/FragmentationAction/StoreSaturatedFragmentAction.cpp

-              r2a03b0
+              ree4f2d
 #include "CodePatterns/Log.hpp"
 #include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp"
+#include "Fragmentation/Exporters/SaturatedFragment.hpp"
 #include "Fragmentation/Graph.hpp"
 #include "World.hpp"
 …
   // store molecule's fragment to file
+  {
+    // we use an empty map here such that saturation is done locally
+    SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
     const enum HydrogenSaturation saturation =  params.DoSaturation.get() ? DoSaturate : DontSaturate;
     ExportGraph_ToFiles exporter(TotalGraph, IncludeHydrogen, saturation);
+    ExportGraph_ToFiles exporter(TotalGraph, IncludeHydrogen, saturation, globalsaturationpositions);
     exporter.setPrefix(params.prefix.get());
     exporter.setOutputTypes(params.types.get());

src/Actions/GlobalListOfActions.hpp

r2a03b0	ree4f2d
21	21	(Redo) \
22	22	(GraphUpdateMolecules) \
	23	(GraphCorrectBondDegree) \
23	24	(GraphCreateAdjacency) \
24	25	(GraphDepthFirstSearch) \

src/Actions/Makefile.am

-              r2a03b0
+              ree4f2d
 GRAPHACTIONSOURCE = \
+  Actions/GraphAction/CorrectBondDegreeAction.cpp \
   Actions/GraphAction/CreateAdjacencyAction.cpp \
   Actions/GraphAction/DepthFirstSearchAction.cpp \
 …
   Actions/GraphAction/UpdateMoleculesAction.cpp
 GRAPHACTIONHEADER = \
+  Actions/GraphAction/CorrectBondDegreeAction.hpp \
   Actions/GraphAction/CreateAdjacencyAction.hpp \
   Actions/GraphAction/DepthFirstSearchAction.hpp \
 …
   Actions/GraphAction/UpdateMoleculesAction.hpp
 GRAPHACTIONDEFS = \
+  Actions/GraphAction/CorrectBondDegreeAction.def \
   Actions/GraphAction/CreateAdjacencyAction.def \
   Actions/GraphAction/DepthFirstSearchAction.def \

src/Atom/atom_bondedparticleinfo.hpp

-              r2a03b0
+              ree4f2d
 #include <list>
+#include <vector>
 /****************************************** forward declarations *****************************/
 …
 class BondedParticle;
+#define BondList list<bond::ptr >
+typedef std::list<bond::ptr > BondList;
 /********************************************** declarations *******************************/

src/Element/elements_db.cpp

-              r2a03b0
+              ree4f2d
 const char *HbondangleDB =\
 "# atomicnumber angles for single, double and triple bond (-1 no angle)\n\
        180     -1      -1\n\
        180     131.0   109.2\n\
        180     120     109.47\n\
        180     110     106.67\n\
        180     104.5   -1\n\
       180     120     109.47\n\
       180     -1      -1\n\
       180     -1      -1\n\
       180 -1 -1\n\
       180     120     109.47\n\
       180 -1 -1\n\
       180 -1 -1\n\
+       0       -1      -1\n\
+       0       131.0   109.2\n\
+       0       120     109.47\n\
+       0       110     106.67\n\
+       0       104.5   -1\n\
+      0       120     109.47\n\
+      0       -1      -1\n\
+      0       -1      -1\n\
+      0 -1 -1\n\
+      0       120     109.47\n\
+      0 -1 -1\n\
+      0 -1 -1\n\
 ";

src/Fragmentation/Exporters/ExportGraph.cpp

-              r2a03b0
+              ree4f2d
   const Graph &_graph,
   const enum HydrogenTreatment _treatment,
+  const enum HydrogenSaturation _saturation) :
+  const enum HydrogenSaturation _saturation,
+  const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
   TotalGraph(_graph),
   BondFragments(World::getPointer()),
   treatment(_treatment),
   saturation(_saturation),
+  globalsaturationpositions(_globalsaturationpositions),
   CurrentKeySet(TotalGraph.begin())
+{
 …
           hydrogens,
           treatment,
+          saturation)
+          saturation,
+          globalsaturationpositions)
   );
   // and return
 …
+}
 /** Internal helper to create from each keyset a molecule
+ *
  */
 void ExportGraph::prepareMolecule()
+{
   size_t count = 0;
   for(Graph::const_iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
     KeySet test = (*runner).first;
     LOG(2, "DEBUG: Fragment No." << (*runner).second.first << " with TEFactor "
         << (*runner).second.second << ".");
     BondFragments.insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
     ++count;
+  }
   LOG(1, "INFO: " << count << "/" << BondFragments.ListOfMolecules.size()
       << " fragments generated from the keysets.");
+}
 /** Stores a fragment from \a KeySet into \a molecule.
  * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
  * molecule and adds missing hydrogen where bonds were cut.
  * \param &Leaflet pointer to KeySet structure
  * \param IsAngstroem whether we have Ansgtroem or bohrradius
  * \return pointer to constructed molecule
  */
 molecule * ExportGraph::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
+{
   Info info(__func__);
   ListOfLocalAtoms_t SonList;
   molecule *Leaf = World::getInstance().createMolecule();
   StoreFragmentFromKeySet_Init(Leaf, Leaflet, SonList);
   // create the bonds between all: Make it an induced subgraph and add hydrogen
 //  LOG(2, "Creating bonds from father graph (i.e. induced subgraph creation).");
   CreateInducedSubgraphOfFragment(Leaf, SonList, IsAngstroem);
   //Leaflet->Leaf->ScanForPeriodicCorrection(out);
   return Leaf;
+}
 /** Initializes some value for putting fragment of \a *mol into \a *Leaf.
  * \param *Leaf fragment molecule
  * \param &Leaflet pointer to KeySet structure
  * \param SonList calloc'd list which atom of \a *Leaf is a son of which atom in \a *mol
  * \return number of atoms in fragment
  */
 int ExportGraph::StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList)
+{
   atom *FatherOfRunner = NULL;
   // first create the minimal set of atoms from the KeySet
   World &world = World::getInstance();
   int size = 0;
   for(KeySet::const_iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
     FatherOfRunner = world.getAtom(AtomById(*runner));  // find the id
     SonList.insert( std::make_pair(FatherOfRunner->getNr(), Leaf->AddCopyAtom(FatherOfRunner) ) );
     size++;
+  }
   return size;
+}
 /** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
  * \param *Leaf fragment molecule
  * \param IsAngstroem whether we have Ansgtroem or bohrradius
  * \param SonList list which atom of \a *Leaf is another atom's son
  */
 void ExportGraph::CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem)
+{
   bool LonelyFlag = false;
   atom *OtherFather = NULL;
   atom *FatherOfRunner = NULL;
   // we increment the iter just before skipping the hydrogen
   // as we use AddBond, we cannot have a const_iterator here
   for (molecule::iterator iter = Leaf->begin(); iter != Leaf->end();) {
     LonelyFlag = true;
     FatherOfRunner = (*iter)->father;
     ASSERT(FatherOfRunner,"Atom without father found");
     if (SonList.find(FatherOfRunner->getNr()) != SonList.end())  {  // check if this, our father, is present in list
       // create all bonds
       const BondList& ListOfBonds = FatherOfRunner->getListOfBonds();
       for (BondList::const_iterator BondRunner = ListOfBonds.begin();
           BondRunner != ListOfBonds.end();
           ++BondRunner) {
         OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
         if (SonList.find(OtherFather->getNr()) != SonList.end()) {
 //          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
 //              << " is bound to " << *OtherFather << ", whose son is "
 //              << *SonList[OtherFather->getNr()] << ".");
           if (OtherFather->getNr() > FatherOfRunner->getNr()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
             std::stringstream output;
 //            output << "ACCEPT: Adding Bond: "
             output << Leaf->AddBond((*iter), SonList[OtherFather->getNr()], (*BondRunner)->getDegree());
 //            LOG(3, output.str());
             //NumBonds[(*iter)->getNr()]++;
           } else {
 //            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
+          }
           LonelyFlag = false;
         } else {
 //          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
 //              << " is bound to " << *OtherFather << ", who has no son in this fragment molecule.");
           if (saturation == DoSaturate) {
 //          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
             if (!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
               exit(1);
           } else if ((treatment == ExcludeHydrogen) && (OtherFather->getElementNo() == (atomicNumber_t)1)) {
             // just copy the atom if it's a hydrogen
             atom * const OtherWalker = Leaf->AddCopyAtom(OtherFather);
             Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
+          }
           //NumBonds[(*iter)->getNr()] += Binder->getDegree();
+        }
+      }
     } else {
       ELOG(1, "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->getNr()] << "!");
+    }
     if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
       LOG(0, **iter << "has got bonds only to hydrogens!");
+    }
     ++iter;
     if (saturation == DoSaturate) {
       while ((iter != Leaf->end()) && ((*iter)->getType()->getAtomicNumber() == 1)){ // skip added hydrogen
         iter++;
+      }
+    }
+  }
+}
+///** Internal helper to create from each keyset a molecule
+// *
+// */
+//void ExportGraph::prepareMolecule()
+//{
+//  size_t count = 0;
+//  for(Graph::const_iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
+//    KeySet test = (*runner).first;
+//    LOG(2, "DEBUG: Fragment No." << (*runner).second.first << " with TEFactor "
+//        << (*runner).second.second << ".");
+//    BondFragments.insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
+//    ++count;
+//  }
+//  LOG(1, "INFO: " << count << "/" << BondFragments.ListOfMolecules.size()
+//      << " fragments generated from the keysets.");
+//}
+//
+///** Stores a fragment from \a KeySet into \a molecule.
+// * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
+// * molecule and adds missing hydrogen where bonds were cut.
+// * \param &Leaflet pointer to KeySet structure
+// * \param IsAngstroem whether we have Ansgtroem or bohrradius
+// * \return pointer to constructed molecule
+// */
+//molecule * ExportGraph::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
+//{
+//  Info info(__func__);
+//  ListOfLocalAtoms_t SonList;
+//  molecule *Leaf = World::getInstance().createMolecule();
+//
+//  StoreFragmentFromKeySet_Init(Leaf, Leaflet, SonList);
+//  // create the bonds between all: Make it an induced subgraph and add hydrogen
+////  LOG(2, "Creating bonds from father graph (i.e. induced subgraph creation).");
+//  CreateInducedSubgraphOfFragment(Leaf, SonList, IsAngstroem);
+//
+//  //Leaflet->Leaf->ScanForPeriodicCorrection(out);
+//  return Leaf;
+//}
+//
+///** Initializes some value for putting fragment of \a *mol into \a *Leaf.
+// * \param *Leaf fragment molecule
+// * \param &Leaflet pointer to KeySet structure
+// * \param SonList calloc'd list which atom of \a *Leaf is a son of which atom in \a *mol
+// * \return number of atoms in fragment
+// */
+//int ExportGraph::StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList)
+//{
+//  atom *FatherOfRunner = NULL;
+//
+//  // first create the minimal set of atoms from the KeySet
+//  World &world = World::getInstance();
+//  int size = 0;
+//  for(KeySet::const_iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
+//    FatherOfRunner = world.getAtom(AtomById(*runner));  // find the id
+//    SonList.insert( std::make_pair(FatherOfRunner->getNr(), Leaf->AddCopyAtom(FatherOfRunner) ) );
+//    size++;
+//  }
+//  return size;
+//}
+//
+///** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
+// * \param *Leaf fragment molecule
+// * \param IsAngstroem whether we have Ansgtroem or bohrradius
+// * \param SonList list which atom of \a *Leaf is another atom's son
+// */
+//void ExportGraph::CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem)
+//{
+//  bool LonelyFlag = false;
+//  atom *OtherFather = NULL;
+//  atom *FatherOfRunner = NULL;
+//
+//  // we increment the iter just before skipping the hydrogen
+//  // as we use AddBond, we cannot have a const_iterator here
+//  for (molecule::iterator iter = Leaf->begin(); iter != Leaf->end();) {
+//    LonelyFlag = true;
+//    FatherOfRunner = (*iter)->father;
+//    ASSERT(FatherOfRunner,"Atom without father found");
+//    if (SonList.find(FatherOfRunner->getNr()) != SonList.end())  {  // check if this, our father, is present in list
+//      // create all bonds
+//      const BondList& ListOfBonds = FatherOfRunner->getListOfBonds();
+//      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+//          BondRunner != ListOfBonds.end();
+//          ++BondRunner) {
+//        OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
+//        if (SonList.find(OtherFather->getNr()) != SonList.end()) {
+////          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
+////              << " is bound to " << *OtherFather << ", whose son is "
+////              << *SonList[OtherFather->getNr()] << ".");
+//          if (OtherFather->getNr() > FatherOfRunner->getNr()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
+//            std::stringstream output;
+////            output << "ACCEPT: Adding Bond: "
+//            output << Leaf->AddBond((*iter), SonList[OtherFather->getNr()], (*BondRunner)->getDegree());
+////            LOG(3, output.str());
+//            //NumBonds[(*iter)->getNr()]++;
+//          } else {
+////            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
+//          }
+//          LonelyFlag = false;
+//        } else {
+////          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
+////              << " is bound to " << *OtherFather << ", who has no son in this fragment molecule.");
+//          if (saturation == DoSaturate) {
+////          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
+//            if (!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
+//              exit(1);
+//          } else if ((treatment == ExcludeHydrogen) && (OtherFather->getElementNo() == (atomicNumber_t)1)) {
+//            // just copy the atom if it's a hydrogen
+//            atom * const OtherWalker = Leaf->AddCopyAtom(OtherFather);
+//            Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
+//          }
+//          //NumBonds[(*iter)->getNr()] += Binder->getDegree();
+//        }
+//      }
+//    } else {
+//      ELOG(1, "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->getNr()] << "!");
+//    }
+//    if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
+//      LOG(0, **iter << "has got bonds only to hydrogens!");
+//    }
+//    ++iter;
+//    if (saturation == DoSaturate) {
+//      while ((iter != Leaf->end()) && ((*iter)->getType()->getAtomicNumber() == 1)){ // skip added hydrogen
+//        iter++;
+//      }
+//    }
+//  }
+//}

src/Fragmentation/Exporters/ExportGraph.hpp

-              r2a03b0
+              ree4f2d
       const Graph &_graph,
       const enum HydrogenTreatment _treatment,
+      const enum HydrogenSaturation _saturation);
+      const enum HydrogenSaturation _saturation,
+      const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
   virtual ~ExportGraph();
 …
   void releaseFragment(SaturatedFragment_ptr &_ptr);
   void prepareMolecule();
   molecule * StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem);
   int StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList);
   void CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem);
+//  void prepareMolecule();
+//  molecule * StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem);
+//  int StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList);
+//  void CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem);
 protected:
 …
   const enum HydrogenSaturation saturation;
+  //!> Global information over all atoms with saturation positions to be used per fragment
+  const SaturatedFragment::GlobalSaturationPositions_t &globalsaturationpositions;
 private:
   //!> iterator pointing at the CurrentKeySet to be exported

src/Fragmentation/Exporters/ExportGraph_ToFiles.cpp

-              r2a03b0
+              ree4f2d
  * @param _treatment whether to always add already present hydrogens or not
  * @param _saturation whether to saturate dangling bonds with hydrogen or not
+ * @param _globalsaturationpositions possibly empty map with global information
+ *        where to place saturation hydrogens to fulfill consistency principle
  */
 ExportGraph_ToFiles::ExportGraph_ToFiles(
     const Graph &_graph,
     const enum HydrogenTreatment _treatment,
+    const enum HydrogenSaturation _saturation) :
+                ExportGraph(_graph, _treatment, _saturation)
+    const enum HydrogenSaturation _saturation,
+    const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
+                ExportGraph(_graph, _treatment, _saturation, _globalsaturationpositions)
 {}

src/Fragmentation/Exporters/ExportGraph_ToFiles.hpp

-              r2a03b0
+              ree4f2d
             const Graph &_graph,
             const enum HydrogenTreatment _treatment,
+            const enum HydrogenSaturation _saturation);
+            const enum HydrogenSaturation _saturation,
+      const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
         virtual ~ExportGraph_ToFiles();

src/Fragmentation/Exporters/ExportGraph_ToJobs.cpp

-              r2a03b0
+              ree4f2d
     const Graph &_graph,
     const enum HydrogenTreatment _treatment,
+    const enum HydrogenSaturation _saturation) :
+    ExportGraph(_graph, _treatment, _saturation),
+    const enum HydrogenSaturation _saturation,
+    const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
+    ExportGraph(_graph, _treatment, _saturation,_globalsaturationpositions),
     level(5)
 {}

src/Fragmentation/Exporters/ExportGraph_ToJobs.hpp

-              r2a03b0
+              ree4f2d
    * \param _treatment whether hydrogen is excluded in the _graph or not
    * \param _saturation whether we saturate dangling bonds or not
+   * \param _globalsaturationpositions possibly empty map with global information
+   *        where to place saturation hydrogens to fulfill consistency principle
    */
   ExportGraph_ToJobs(
       const Graph &_graph,
       const enum HydrogenTreatment _treatment,
+      const enum HydrogenSaturation _saturation);
+      const enum HydrogenSaturation _saturation,
+      const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
   virtual ~ExportGraph_ToJobs();

src/Fragmentation/Exporters/SaturatedFragment.cpp

-              r2a03b0
+              ree4f2d
     HydrogenPool &_hydrogens,
     const enum HydrogenTreatment _treatment,
+    const enum HydrogenSaturation _saturation) :
+    const enum HydrogenSaturation _saturation,
+    const GlobalSaturationPositions_t &_globalsaturationpositions) :
   container(_container),
   set(_set),
 …
   container.insert(set);
+  // prepare saturation hydrogens
+  saturate();
+  // prepare saturation hydrogens, either using global information
+  // or if not given, local information (created in the function)
+  if (_globalsaturationpositions.empty())
+    saturate();
+  else
+    saturate(_globalsaturationpositions);
+}
 …
+}
+void SaturatedFragment::saturate()
+{
+  // gather all atoms in a vector
+  std::vector<atom *> atoms;
+  for (KeySet::const_iterator iter = FullMolecule.begin();
+      iter != FullMolecule.end();
+typedef std::vector<atom *> atoms_t;
+atoms_t gatherAllAtoms(const KeySet &_FullMolecule)
+{
+  atoms_t atoms;
+  for (KeySet::const_iterator iter = _FullMolecule.begin();
+      iter != _FullMolecule.end();
       ++iter) {
     atom * const Walker = World::getInstance().getAtom(AtomById(*iter));
     ASSERT( Walker != NULL,
         "SaturatedFragment::OutputConfig() - id "
+        "gatherAllAtoms() - id "
         +toString(*iter)+" is unknown to World.");
     atoms.push_back(Walker);
+  }
+//  bool LonelyFlag = false;
+  for (std::vector<atom *>::const_iterator iter = atoms.begin();
+      iter != atoms.end();
+  return atoms;
+}
+typedef std::map<atom *, BondList > CutBonds_t;
+CutBonds_t gatherCutBonds(
+    const atoms_t &_atoms,
+    const KeySet &_set,
+    const enum HydrogenTreatment _treatment)
+{
+  //  bool LonelyFlag = false;
+  CutBonds_t CutBonds;
+  for (atoms_t::const_iterator iter = _atoms.begin();
+      iter != _atoms.end();
       ++iter) {
     atom * const Walker = *iter;
 …
         ++BondRunner) {
       atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
       // if in set
       if (set.find(OtherWalker->getId()) != set.end()) {
+      // if other atom is in key set or is a specially treated hydrogen
+      if (_set.find(OtherWalker->getId()) != _set.end()) {
         LOG(4, "DEBUG: Walker " << *Walker << " is bound to " << *OtherWalker << ".");
+//        if (OtherWalker->getId() > Walker->getId()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
+////          std::stringstream output;
+////            output << "ACCEPT: Adding Bond: "
+//          output << Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
+////            LOG(3, output.str());
+//          //NumBonds[(*iter)->getNr()]++;
+//        } else {
+////            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
+//        }
+//        LonelyFlag = false;
+      } else if ((_treatment == ExcludeHydrogen)
+          && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
+        LOG(4, "DEBUG: Walker " << *Walker << " is bound to specially treated hydrogen " <<
+            *OtherWalker << ".");
       } else {
         LOG(4, "DEBUG: Walker " << *Walker << " is bound to "
             << *OtherWalker << ", who is not in this fragment molecule.");
+        if (saturation == DoSaturate) {
+//          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
+          if (!AddHydrogenReplacementAtom(
+              (*BondRunner),
+              Walker,
+              OtherWalker,
+              World::getInstance().getConfig()->IsAngstroem == 1))
+            exit(1);
+          if (CutBonds.count(Walker) == 0)
+            CutBonds.insert( std::make_pair(Walker, BondList() ));
+          CutBonds[Walker].push_back(*BondRunner);
+      }
+    }
+  }
+  return CutBonds;
+}
+typedef std::vector<atomId_t> atomids_t;
+atomids_t gatherPresentExcludedHydrogens(
+    const atoms_t &_atoms,
+    const KeySet &_set,
+    const enum HydrogenTreatment _treatment)
+{
+  //  bool LonelyFlag = false;
+  atomids_t ExcludedHydrogens;
+  for (atoms_t::const_iterator iter = _atoms.begin();
+      iter != _atoms.end();
+      ++iter) {
+    atom * const Walker = *iter;
+    // go through all bonds
+    const BondList& ListOfBonds = Walker->getListOfBonds();
+    for (BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner) {
+      atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
+      // if other atom is in key set or is a specially treated hydrogen
+      if (_set.find(OtherWalker->getId()) != _set.end()) {
+        LOG(6, "DEBUG: OtherWalker " << *OtherWalker << " is in set already.");
+      } else if ((_treatment == ExcludeHydrogen)
+          && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
+        LOG(5, "DEBUG: Adding excluded hydrogen OtherWalker " << *OtherWalker << ".");
+        ExcludedHydrogens.push_back(OtherWalker->getId());
+      } else {
+        LOG(6, "DEBUG: OtherWalker " << *Walker << " is not in this fragment molecule and no hydrogen.");
+      }
+    }
+  }
+  return ExcludedHydrogens;
+}
+void SaturatedFragment::saturate()
+{
+  // so far, we just have a set of keys. Hence, convert to atom refs
+  // and gather all atoms in a vector
+  std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
+  // go through each atom of the fragment and gather all cut bonds in list
+  CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
+  // add excluded hydrogens to FullMolecule if treated specially
+  if (treatment == ExcludeHydrogen) {
+    atomids_t ExcludedHydrogens = gatherPresentExcludedHydrogens(atoms, set, treatment);
+    FullMolecule.insert(ExcludedHydrogens.begin(), ExcludedHydrogens.end());
+  }
+  // go through all cut bonds and replace with a hydrogen
+  for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
+      atomiter != CutBonds.end(); ++atomiter) {
+    atom * const Walker = atomiter->first;
+    if (!saturateAtom(Walker, atomiter->second))
+      exit(1);
+  }
+}
+void SaturatedFragment::saturate(
+    const GlobalSaturationPositions_t &_globalsaturationpositions)
+{
+  // so far, we just have a set of keys. Hence, convert to atom refs
+  // and gather all atoms in a vector
+  std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
+  // go through each atom of the fragment and gather all cut bonds in list
+  CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
+  // add excluded hydrogens to FullMolecule if treated specially
+  if (treatment == ExcludeHydrogen) {
+    atomids_t ExcludedHydrogens = gatherPresentExcludedHydrogens(atoms, set, treatment);
+    FullMolecule.insert(ExcludedHydrogens.begin(), ExcludedHydrogens.end());
+  }
+  // go through all cut bonds and replace with a hydrogen
+  if (saturation == DoSaturate) {
+    for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
+        atomiter != CutBonds.end(); ++atomiter) {
+      atom * const Walker = atomiter->first;
+      LOG(4, "DEBUG: We are now saturating dangling bonds of " << *Walker);
+      // gather set of positions for this atom from global map
+      GlobalSaturationPositions_t::const_iterator mapiter =
+          _globalsaturationpositions.find(Walker->getId());
+      ASSERT( mapiter != _globalsaturationpositions.end(),
+          "SaturatedFragment::saturate() - no global information for "
+          +toString(*Walker));
+      const SaturationsPositionsPerNeighbor_t &saturationpositions =
+          mapiter->second;
+      // go through all cut bonds for this atom
+      for (BondList::const_iterator bonditer = atomiter->second.begin();
+          bonditer != atomiter->second.end(); ++bonditer) {
+        atom * const OtherWalker = (*bonditer)->GetOtherAtom(Walker);
+        // get positions from global map
+        SaturationsPositionsPerNeighbor_t::const_iterator positionsiter =
+            saturationpositions.find(OtherWalker->getId());
+        ASSERT(positionsiter != saturationpositions.end(),
+            "SaturatedFragment::saturate() - no information on bond neighbor "
+            +toString(*OtherWalker)+" to atom "+toString(*Walker));
+        ASSERT(!positionsiter->second.empty(),
+            "SaturatedFragment::saturate() - no positions for saturating bond to"
+            +toString(*OtherWalker)+" to atom "+toString(*Walker));
+//        // get typical bond distance from elements database
+//        double BondDistance = Walker->getType()->getHBondDistance(positionsiter->second.size()-1);
+//        if (BondDistance < 0.) {
+//          ELOG(2, "saturateAtoms() - no typical hydrogen bond distance of degree "
+//              +toString(positionsiter->second.size())+" for element "
+//              +toString(Walker->getType()->getName()));
+//          // try bond degree 1 distance
+//          BondDistance = Walker->getType()->getHBondDistance(1-1);
+//          if (BondDistance < 0.) {
+//            ELOG(1, "saturateAtoms() - no typical hydrogen bond distance for element "
+//                +toString(Walker->getType()->getName()));
+//            BondDistance = 1.;
+//          }
+//        }
+        // place hydrogen at each point
+        LOG(4, "DEBUG: Places to saturate for atom " << *OtherWalker
+            << " are " << positionsiter->second);
+        atom * const father = Walker;
+        for (SaturationsPositions_t::const_iterator positer = positionsiter->second.begin();
+            positer != positionsiter->second.end(); ++positer) {
+          const atom& hydrogen =
+              setHydrogenReplacement(Walker, *positer, 1., father);
+          FullMolecule.insert(hydrogen.getId());
+        }
-//        } else if ((treatment == ExcludeHydrogen) && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
-//          // just copy the atom if it's a hydrogen
-//          atom * const OtherWalker = Leaf->AddCopyAtom(OtherWalker);
-//          Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
-//        }
-        //NumBonds[(*iter)->getNr()] += Binder->getDegree();
+      }
+    }
+  }
+  } else
+    LOG(3, "INFO: We are not saturating cut bonds.");
+}
+const atom& SaturatedFragment::setHydrogenReplacement(
+  const atom * const _OwnerAtom,
+  const Vector &_position,
+  const double _distance,
+  atom * const _father)
+{
+  atom * const _atom = hydrogens.leaseHydrogen();    // new atom
+  _atom->setPosition( _OwnerAtom->getPosition() + _distance * _position );
+  // always set as fixed ion (not moving during molecular dynamics simulation)
+  _atom->setFixedIon(true);
+  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
+  _atom->father = _father;
+  SaturationHydrogens.insert(_atom->getId());
+  return *_atom;
+}
+bool SaturatedFragment::saturateAtom(
+    atom * const _atom,
+    const BondList &_cutbonds)
+{
+  // go through each bond and replace
+  for (BondList::const_iterator bonditer = _cutbonds.begin();
+      bonditer != _cutbonds.end(); ++bonditer) {
+    atom * const OtherWalker = (*bonditer)->GetOtherAtom(_atom);
+    if (!AddHydrogenReplacementAtom(
+        (*bonditer),
+        _atom,
+        OtherWalker,
+        World::getInstance().getConfig()->IsAngstroem == 1))
+      return false;
+  }
+  return true;
+}
 …
   if (BondRescale == -1) {
     ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->getDegree() << "!");
+    return false;
+    BondRescale = bondlength;
+    BondRescale = Origin->getType()->getHBondDistance(TopBond->getDegree());
+    if (BondRescale == -1) {
+      ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of any degree!");
+      return false;
+      BondRescale = bondlength;
+    }
   } else {
     if (!IsAngstroem)

src/Fragmentation/Exporters/SaturatedFragment.hpp

-              r2a03b0
+              ree4f2d
 #include <string>
+#include "Atom/atom_bondedparticleinfo.hpp"
 #include "Bond/bond.hpp"
 #include "Fragmentation/KeySet.hpp"
 #include "Fragmentation/HydrogenSaturation_enum.hpp"
 #include "Parser/FormatParserStorage.hpp"
+#include "LinearAlgebra/Vector.hpp"
 class atom;
 …
   typedef std::set<KeySet> KeySetsInUse_t;
+  //!> List of points giving saturation positions for a single bond neighbor
+  typedef std::list<Vector> SaturationsPositions_t;
+  //!> map for one atom, containing the saturation points for all its neighbors
+  typedef std::map<int, SaturationsPositions_t> SaturationsPositionsPerNeighbor_t;
+  //!> containing the saturation points over all desired atoms required
+  typedef std::map<int, SaturationsPositionsPerNeighbor_t> GlobalSaturationPositions_t;
   /** Constructor of SaturatedFragment requires \a set which we are tightly
    * associated.
 …
    * \param _container container to add KeySet as in-use
    * \param _hydrogens pool with hydrogens for saturation
+   * \param _globalsaturationpositions saturation positions to be used
    */
   SaturatedFragment(
 …
       HydrogenPool &_hydrogens,
       const enum HydrogenTreatment _treatment,
+      const enum HydrogenSaturation saturation);
+      const enum HydrogenSaturation saturation,
+      const GlobalSaturationPositions_t &_globalsaturationpositions);
   /** Destructor of class SaturatedFragment.
 …
   /** Helper function to lease and bring in place saturation hydrogens.
+   *
+   * Here, we use local information to calculate saturation positions.
+   *
    */
   void saturate();
+  /** Helper function to lease and bring in place saturation hydrogens.
+   *
+   * Here, saturation positions have to be calculated before and are fully
+   * stored in \a _globalsaturationpositions.
+   *
+   * \param_globalsaturationpositions
+   */
+  void saturate(const GlobalSaturationPositions_t &_globalsaturationpositions);
+  /** Replaces all cut bonds with respect to the given atom by hydrogens.
+   *
+   * \param _atom atom whose cut bonds to saturate
+   * \param _cutbonds list of cut bonds for \a _atom
+   * \return true - bonds saturated, false - something went wrong
+   */
+  bool saturateAtom(atom * const _atom, const BondList &_cutbonds);
   /** Helper function to get a hydrogen replacement for a given \a replacement.
 …
    */
   atom * const getHydrogenReplacement(atom * const replacement);
+  /** Sets a saturation hydrogen at the given position in place of \a _father.
+   *
+   * \param _OwnerAtom atom "owning" the hydrogen (i.e. the one getting saturated)
+   * \param _position new position relative to \a _OwnerAtom
+   * \param _distance scale factor to the distance (default 1.)
+   * \param _father bond partner of \a _OwnerAtom that is replaced
+   */
+  const atom& setHydrogenReplacement(
+      const atom * const _OwnerAtom,
+      const Vector &_position,
+      const double _distance,
+      atom * const _father);
   /** Leases and adds a Hydrogen atom in replacement for the given atom \a *partner in bond with a *origin.

src/Fragmentation/Exporters/unittests/Makefile.am

-              r2a03b0
+              ree4f2d
 FRAGMENTATIONEXPORTERSSOURCES = \
         ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.cpp \
+        ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp
+        ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp \
+        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.cpp
 FRAGMENTATIONEXPORTERSTESTSHEADERS = \
         ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.hpp \
+        ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.hpp
+        ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.hpp \
+        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.hpp
 FRAGMENTATIONEXPORTERSTESTS = \
         HydrogenPoolUnitTest \
+        SaturatedFragmentUnitTest
+        SaturatedFragmentUnitTest \
+        SaturationDistanceMaximizerUnitTest
 TESTS += $(FRAGMENTATIONEXPORTERSTESTS)
 …
         $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
         ${FRAGMENTATIONEXPORTERSLIBS}
+SaturationDistanceMaximizerUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.cpp \
+        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.hpp \
+        ../Fragmentation/Exporters/unittests/stubs/SaturatedBondStub.cpp
+SaturationDistanceMaximizerUnitTest_LDADD = \
+        ../libMolecuilderUI.la  \
+        $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
+        ${FRAGMENTATIONEXPORTERSLIBS}

src/Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp

-              r2a03b0
+              ree4f2d
   ASSERT_DO(Assert::Throw);
+  SaturatedFragment::GlobalSaturationPositions_t globalpositions;
   set = new KeySet;
+  fragment = new SaturatedFragment(*set, KeySetsInUse, hydrogens, ExcludeHydrogen, DoSaturate);
+  fragment = new SaturatedFragment(
+      *set,
+      KeySetsInUse,
+      hydrogens,
+      ExcludeHydrogen,
+      DoSaturate,
+      globalpositions);
+}

src/Fragmentation/Makefile.am

-              r2a03b0
+              ree4f2d
         Fragmentation/Exporters/ExportGraph.cpp \
         Fragmentation/Exporters/HydrogenPool.cpp \
+        Fragmentation/Exporters/SaturatedBond.cpp \
         Fragmentation/Exporters/SaturatedFragment.cpp \
+        Fragmentation/Exporters/SaturationDistanceMaximizer.cpp \
         Fragmentation/Homology/FragmentEdge.cpp \
         Fragmentation/Homology/FragmentNode.cpp \
 …
         Fragmentation/Exporters/ExportGraph.hpp \
         Fragmentation/Exporters/HydrogenPool.hpp \
+        Fragmentation/Exporters/SaturatedBond.hpp \
         Fragmentation/Exporters/SaturatedFragment.hpp \
+        Fragmentation/Exporters/SaturationDistanceMaximizer.hpp \
         Fragmentation/Homology/FragmentEdge.hpp \
         Fragmentation/Homology/FragmentNode.hpp \

src/Fragmentation/Summation/Converter/DataConverter.hpp

-              r2a03b0
+              ree4f2d
     MPQCDataForceMap_t instance;
     // must convert int to index_t
+    if (DoLog(5)) {
+      std::stringstream output;
+      for (KeySetsContainer::IntVector::const_iterator outiter = arrayiter->begin();
+          outiter != arrayiter->end(); ++outiter) {
+        output << *outiter << "\t";
+      }
+      LOG(5, "DEBUG: indices are " << output.str());
+    }
     IndexedVectors::indices_t indices(arrayiter->begin(), arrayiter->end());
     boost::fusion::at_key<MPQCDataFused::forces>(instance) =

src/Graph/Makefile.am

-              r2a03b0
+              ree4f2d
 GRAPHSOURCE = \
         Graph/BondGraph.cpp \
-        Graph/BreadthFirstSearchAdd.cpp \
         Graph/BuildInducedSubgraph.cpp \
         Graph/AdjacencyList.cpp \
 …
 GRAPHHEADER = \
         Graph/BondGraph.hpp \
-        Graph/BreadthFirstSearchAdd.hpp \
         Graph/BuildInducedSubgraph.hpp \
         Graph/AdjacencyList.hpp \

src/Makefile.am

-              r2a03b0
+              ree4f2d
 include RandomNumbers/Makefile.am
 include Shapes/Makefile.am
+include Tesselation/Makefile.am
 include UIElements/Makefile.am
 …
   Thermostats/Woodcock.hpp
-TESSELATIONSOURCE = \
-  Tesselation/ApproximateShapeArea.cpp \
-  Tesselation/ApproximateShapeVolume.cpp \
-  Tesselation/boundary.cpp \
-  Tesselation/BoundaryLineSet.cpp \
-  Tesselation/BoundaryPointSet.cpp \
-  Tesselation/BoundaryPolygonSet.cpp \
-  Tesselation/BoundaryTriangleSet.cpp \
-  Tesselation/CandidateForTesselation.cpp \
-  Tesselation/ellipsoid.cpp \
-  Tesselation/tesselation.cpp \
-  Tesselation/tesselationhelpers.cpp \
-  Tesselation/triangleintersectionlist.cpp
-TESSELATIONHEADER = \
-  Tesselation/ApproximateShapeArea.hpp \
-  Tesselation/ApproximateShapeVolume.hpp \
-  Tesselation/boundary.hpp \
-  Tesselation/BoundaryLineSet.hpp \
-  Tesselation/BoundaryMaps.hpp \
-  Tesselation/BoundaryPointSet.hpp \
-  Tesselation/BoundaryPolygonSet.hpp \
-  Tesselation/BoundaryTriangleSet.hpp \
-  Tesselation/CandidateForTesselation.hpp \
-  Tesselation/ellipsoid.hpp \
-  Tesselation/tesselation.hpp \
-  Tesselation/tesselationhelpers.hpp \
-  Tesselation/triangleintersectionlist.hpp
 MOLECUILDERSOURCE = \
   ${BONDSOURCE} \
 …
   ${DYNAMICSSOURCE} \
   ${THERMOSTATSOURCE} \
-  ${TESSELATIONSOURCE} \
   Shapes/ShapeFactory.cpp \
   AtomIdSet.cpp \
 …
   ${DYNAMICSHEADER} \
   ${THERMOSTATHEADER} \
-  ${TESSELATIONHEADER} \
   Shapes/ShapeFactory.hpp \
   AtomIdSet.hpp \
 …
         $(BOOST_THREAD_LDFLAGS)
 libMolecuilder_la_LIBADD = \
+        libMolecuilderTesselation.la \
         libMolecuilderShapes.la \
         $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \

src/molecule.cpp

-              r2a03b0
+              ree4f2d
  * \todo double and triple bonds splitting (always use the tetraeder angle!)
  */
 bool molecule::AddHydrogenReplacementAtom(bond::ptr TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
+{
 //  Info info(__func__);
   bool AllWentWell = true;    // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
   double bondlength;  // bond length of the bond to be replaced/cut
   double bondangle;  // bond angle of the bond to be replaced/cut
   double BondRescale;   // rescale value for the hydrogen bond length
   bond::ptr FirstBond;
   bond::ptr SecondBond; // Other bonds in double bond case to determine "other" plane
   atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
   double b,l,d,f,g, alpha, factors[NDIM];    // hold temporary values in triple bond case for coordination determination
   Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
   Vector InBondvector;    // vector in direction of *Bond
   const RealSpaceMatrix &matrix =  World::getInstance().getDomain().getM();
   bond::ptr Binder;
   // create vector in direction of bond
   InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
   bondlength = InBondvector.Norm();
    // is greater than typical bond distance? Then we have to correct periodically
    // the problem is not the H being out of the box, but InBondvector have the wrong direction
    // due to TopReplacement or Origin being on the wrong side!
   const BondGraph * const BG = World::getInstance().getBondGraph();
   const range<double> MinMaxBondDistance(
       BG->getMinMaxDistance(TopOrigin,TopReplacement));
   if (!MinMaxBondDistance.isInRange(bondlength)) {
 //    LOG(4, "InBondvector is: " << InBondvector << ".");
     Orthovector1.Zero();
     for (int i=NDIM;i--;) {
       l = TopReplacement->at(i) - TopOrigin->at(i);
       if (fabs(l) > MinMaxBondDistance.last) { // is component greater than bond distance (check against min not useful here)
         Orthovector1[i] = (l < 0) ? -1. : +1.;
       } // (signs are correct, was tested!)
+    }
     Orthovector1 *= matrix;
     InBondvector -= Orthovector1; // subtract just the additional translation
     bondlength = InBondvector.Norm();
 //    LOG(4, "INFO: Corrected InBondvector is now: " << InBondvector << ".");
   } // periodic correction finished
   InBondvector.Normalize();
   // get typical bond length and store as scale factor for later
   ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
   BondRescale = TopOrigin->getType()->getHBondDistance(TopBond->getDegree()-1);
   if (BondRescale == -1) {
     ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
     return false;
     BondRescale = bondlength;
   } else {
     if (!IsAngstroem)
       BondRescale /= (1.*AtomicLengthToAngstroem);
+  }
   // discern single, double and triple bonds
   switch(TopBond->getDegree()) {
     case 1:
       FirstOtherAtom = World::getInstance().createAtom();    // new atom
       FirstOtherAtom->setType(1);  // element is Hydrogen
       FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
       FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
         FirstOtherAtom->father = TopReplacement;
         BondRescale = bondlength;
       } else {
         FirstOtherAtom->father = NULL;  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
+      }
       InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
       FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
 //      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
       Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
       Binder->Cyclic = false;
       Binder->Type = GraphEdge::TreeEdge;
       break;
     case 2:
+      {
         // determine two other bonds (warning if there are more than two other) plus valence sanity check
         const BondList& ListOfBonds = TopOrigin->getListOfBonds();
         for (BondList::const_iterator Runner = ListOfBonds.begin();
             Runner != ListOfBonds.end();
             ++Runner) {
           if ((*Runner) != TopBond) {
             if (FirstBond == NULL) {
               FirstBond = (*Runner);
               FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
             } else if (SecondBond == NULL) {
               SecondBond = (*Runner);
               SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
             } else {
               ELOG(2, "Detected more than four bonds for atom " << TopOrigin->getName());
+            }
+          }
+        }
+      }
       if (SecondOtherAtom == NULL) {  // then we have an atom with valence four, but only 3 bonds: one to replace and one which is TopBond (third is FirstBond)
         SecondBond = TopBond;
         SecondOtherAtom = TopReplacement;
+      }
       if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
 //        LOG(3, "Regarding the double bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << TopOrigin->Name << " to determine orthogonal plane.");
         // determine the plane of these two with the *origin
         try {
           Orthovector1 = Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
+        }
         catch(LinearDependenceException &excp){
           LOG(0, boost::diagnostic_information(excp));
           // TODO: figure out what to do with the Orthovector in this case
           AllWentWell = false;
+        }
       } else {
         Orthovector1.GetOneNormalVector(InBondvector);
+      }
       //LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
       // orthogonal vector and bond vector between origin and replacement form the new plane
       Orthovector1.MakeNormalTo(InBondvector);
       Orthovector1.Normalize();
       //LOG(3, "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << ".");
       // create the two Hydrogens ...
       FirstOtherAtom = World::getInstance().createAtom();
       SecondOtherAtom = World::getInstance().createAtom();
       FirstOtherAtom->setType(1);
       SecondOtherAtom->setType(1);
       FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
       FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
       SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
       SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
       bondangle = TopOrigin->getType()->getHBondAngle(1);
       if (bondangle == -1) {
         ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
         return false;
         bondangle = 0;
+      }
       bondangle *= M_PI/180./2.;
 //      LOG(3, "INFO: ReScaleCheck: InBondvector " << InBondvector << ", " << Orthovector1 << ".");
 //      LOG(3, "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle));
       FirstOtherAtom->Zero();
       SecondOtherAtom->Zero();
       for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
         FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
         SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
+      }
       FirstOtherAtom->Scale(BondRescale);  // rescale by correct BondDistance
       SecondOtherAtom->Scale(BondRescale);
       //LOG(3, "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << ".");
       *FirstOtherAtom += TopOrigin->getPosition();
       *SecondOtherAtom += TopOrigin->getPosition();
       // ... and add to molecule
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
       AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
 //      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
 //      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
       Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
       Binder->Cyclic = false;
       Binder->Type = GraphEdge::TreeEdge;
       Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
       Binder->Cyclic = false;
       Binder->Type = GraphEdge::TreeEdge;
       break;
     case 3:
       // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
       FirstOtherAtom = World::getInstance().createAtom();
       SecondOtherAtom = World::getInstance().createAtom();
       ThirdOtherAtom = World::getInstance().createAtom();
       FirstOtherAtom->setType(1);
       SecondOtherAtom->setType(1);
       ThirdOtherAtom->setType(1);
       FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
       FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
       SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       ThirdOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
       ThirdOtherAtom->setFixedIon(TopReplacement->getFixedIon());
       FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
       SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
       ThirdOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
       // we need to vectors orthonormal the InBondvector
       AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
 //      LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
       try{
         Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
+      }
       catch(LinearDependenceException &excp) {
         LOG(0, boost::diagnostic_information(excp));
         AllWentWell = false;
+      }
 //      LOG(3, "INFO: Orthovector2: " << Orthovector2 << ".")
       // create correct coordination for the three atoms
       alpha = (TopOrigin->getType()->getHBondAngle(2))/180.*M_PI/2.;  // retrieve triple bond angle from database
       l = BondRescale;        // desired bond length
       b = 2.*l*sin(alpha);    // base length of isosceles triangle
       d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.);   // length for InBondvector
       f = b/sqrt(3.);   // length for Orthvector1
       g = b/2.;         // length for Orthvector2
 //      LOG(3, "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", ");
 //      LOG(3, "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", "  << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g));
       factors[0] = d;
       factors[1] = f;
       factors[2] = 0.;
       FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[1] = -0.5*f;
       factors[2] = g;
       SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[2] = -g;
       ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       // rescale each to correct BondDistance
 //      FirstOtherAtom->x.Scale(&BondRescale);
 //      SecondOtherAtom->x.Scale(&BondRescale);
 //      ThirdOtherAtom->x.Scale(&BondRescale);
       // and relative to *origin atom
       *FirstOtherAtom += TopOrigin->getPosition();
       *SecondOtherAtom += TopOrigin->getPosition();
       *ThirdOtherAtom += TopOrigin->getPosition();
       // ... and add to molecule
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
       AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
       AllWentWell = AllWentWell && AddAtom(ThirdOtherAtom);
 //      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
 //      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
 //      LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
       Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
       Binder->Cyclic = false;
       Binder->Type = GraphEdge::TreeEdge;
       Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
       Binder->Cyclic = false;
       Binder->Type = GraphEdge::TreeEdge;
       Binder = AddBond(BottomOrigin, ThirdOtherAtom, 1);
       Binder->Cyclic = false;
       Binder->Type = GraphEdge::TreeEdge;
       break;
     default:
       ELOG(1, "BondDegree does not state single, double or triple bond!");
       AllWentWell = false;
       break;
+  }
   return AllWentWell;
 };
+//bool molecule::AddHydrogenReplacementAtom(bond::ptr TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
+//{
+////  Info info(__func__);
+//  bool AllWentWell = true;    // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
+//  double bondlength;  // bond length of the bond to be replaced/cut
+//  double bondangle;  // bond angle of the bond to be replaced/cut
+//  double BondRescale;   // rescale value for the hydrogen bond length
+//  bond::ptr FirstBond;
+//  bond::ptr SecondBond; // Other bonds in double bond case to determine "other" plane
+//  atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
+//  double b,l,d,f,g, alpha, factors[NDIM];    // hold temporary values in triple bond case for coordination determination
+//  Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
+//  Vector InBondvector;    // vector in direction of *Bond
+//  const RealSpaceMatrix &matrix =  World::getInstance().getDomain().getM();
+//  bond::ptr Binder;
+//
+//  // create vector in direction of bond
+//  InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
+//  bondlength = InBondvector.Norm();
+//
+//   // is greater than typical bond distance? Then we have to correct periodically
+//   // the problem is not the H being out of the box, but InBondvector have the wrong direction
+//   // due to TopReplacement or Origin being on the wrong side!
+//  const BondGraph * const BG = World::getInstance().getBondGraph();
+//  const range<double> MinMaxBondDistance(
+//      BG->getMinMaxDistance(TopOrigin,TopReplacement));
+//  if (!MinMaxBondDistance.isInRange(bondlength)) {
+////    LOG(4, "InBondvector is: " << InBondvector << ".");
+//    Orthovector1.Zero();
+//    for (int i=NDIM;i--;) {
+//      l = TopReplacement->at(i) - TopOrigin->at(i);
+//      if (fabs(l) > MinMaxBondDistance.last) { // is component greater than bond distance (check against min not useful here)
+//        Orthovector1[i] = (l < 0) ? -1. : +1.;
+//      } // (signs are correct, was tested!)
+//    }
+//    Orthovector1 *= matrix;
+//    InBondvector -= Orthovector1; // subtract just the additional translation
+//    bondlength = InBondvector.Norm();
+////    LOG(4, "INFO: Corrected InBondvector is now: " << InBondvector << ".");
+//  } // periodic correction finished
+//
+//  InBondvector.Normalize();
+//  // get typical bond length and store as scale factor for later
+//  ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
+//  BondRescale = TopOrigin->getType()->getHBondDistance(TopBond->getDegree()-1);
+//  if (BondRescale == -1) {
+//    ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
+//    return false;
+//    BondRescale = bondlength;
+//  } else {
+//    if (!IsAngstroem)
+//      BondRescale /= (1.*AtomicLengthToAngstroem);
+//  }
+//
+//  // discern single, double and triple bonds
+//  switch(TopBond->getDegree()) {
+//    case 1:
+//      FirstOtherAtom = World::getInstance().createAtom();    // new atom
+//      FirstOtherAtom->setType(1);  // element is Hydrogen
+//      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
+//      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
+//      if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
+//        FirstOtherAtom->father = TopReplacement;
+//        BondRescale = bondlength;
+//      } else {
+//        FirstOtherAtom->father = NULL;  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
+//      }
+//      InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
+//      FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
+//      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
+////      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
+//      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
+//      Binder->Cyclic = false;
+//      Binder->Type = GraphEdge::TreeEdge;
+//      break;
+//    case 2:
+//      {
+//        // determine two other bonds (warning if there are more than two other) plus valence sanity check
+//        const BondList& ListOfBonds = TopOrigin->getListOfBonds();
+//        for (BondList::const_iterator Runner = ListOfBonds.begin();
+//            Runner != ListOfBonds.end();
+//            ++Runner) {
+//          if ((*Runner) != TopBond) {
+//            if (FirstBond == NULL) {
+//              FirstBond = (*Runner);
+//              FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
+//            } else if (SecondBond == NULL) {
+//              SecondBond = (*Runner);
+//              SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
+//            } else {
+//              ELOG(2, "Detected more than four bonds for atom " << TopOrigin->getName());
+//            }
+//          }
+//        }
+//      }
+//      if (SecondOtherAtom == NULL) {  // then we have an atom with valence four, but only 3 bonds: one to replace and one which is TopBond (third is FirstBond)
+//        SecondBond = TopBond;
+//        SecondOtherAtom = TopReplacement;
+//      }
+//      if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
+////        LOG(3, "Regarding the double bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << TopOrigin->Name << " to determine orthogonal plane.");
+//
+//        // determine the plane of these two with the *origin
+//        try {
+//          Orthovector1 = Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
+//        }
+//        catch(LinearDependenceException &excp){
+//          LOG(0, boost::diagnostic_information(excp));
+//          // TODO: figure out what to do with the Orthovector in this case
+//          AllWentWell = false;
+//        }
+//      } else {
+//        Orthovector1.GetOneNormalVector(InBondvector);
+//      }
+//      //LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
+//      // orthogonal vector and bond vector between origin and replacement form the new plane
+//      Orthovector1.MakeNormalTo(InBondvector);
+//      Orthovector1.Normalize();
+//      //LOG(3, "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << ".");
+//
+//      // create the two Hydrogens ...
+//      FirstOtherAtom = World::getInstance().createAtom();
+//      SecondOtherAtom = World::getInstance().createAtom();
+//      FirstOtherAtom->setType(1);
+//      SecondOtherAtom->setType(1);
+//      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
+//      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
+//      SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
+//      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
+//      FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
+//      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
+//      bondangle = TopOrigin->getType()->getHBondAngle(1);
+//      if (bondangle == -1) {
+//        ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
+//        return false;
+//        bondangle = 0;
+//      }
+//      bondangle *= M_PI/180./2.;
+////      LOG(3, "INFO: ReScaleCheck: InBondvector " << InBondvector << ", " << Orthovector1 << ".");
+////      LOG(3, "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle));
+//      FirstOtherAtom->Zero();
+//      SecondOtherAtom->Zero();
+//      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
+//        FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
+//        SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
+//      }
+//      FirstOtherAtom->Scale(BondRescale);  // rescale by correct BondDistance
+//      SecondOtherAtom->Scale(BondRescale);
+//      //LOG(3, "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << ".");
+//      *FirstOtherAtom += TopOrigin->getPosition();
+//      *SecondOtherAtom += TopOrigin->getPosition();
+//      // ... and add to molecule
+//      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
+//      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
+////      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
+////      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
+//      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
+//      Binder->Cyclic = false;
+//      Binder->Type = GraphEdge::TreeEdge;
+//      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
+//      Binder->Cyclic = false;
+//      Binder->Type = GraphEdge::TreeEdge;
+//      break;
+//    case 3:
+//      // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
+//      FirstOtherAtom = World::getInstance().createAtom();
+//      SecondOtherAtom = World::getInstance().createAtom();
+//      ThirdOtherAtom = World::getInstance().createAtom();
+//      FirstOtherAtom->setType(1);
+//      SecondOtherAtom->setType(1);
+//      ThirdOtherAtom->setType(1);
+//      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
+//      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
+//      SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
+//      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
+//      ThirdOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
+//      ThirdOtherAtom->setFixedIon(TopReplacement->getFixedIon());
+//      FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
+//      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
+//      ThirdOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
+//
+//      // we need to vectors orthonormal the InBondvector
+//      AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
+////      LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
+//      try{
+//        Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
+//      }
+//      catch(LinearDependenceException &excp) {
+//        LOG(0, boost::diagnostic_information(excp));
+//        AllWentWell = false;
+//      }
+////      LOG(3, "INFO: Orthovector2: " << Orthovector2 << ".")
+//
+//      // create correct coordination for the three atoms
+//      alpha = (TopOrigin->getType()->getHBondAngle(2))/180.*M_PI/2.;  // retrieve triple bond angle from database
+//      l = BondRescale;        // desired bond length
+//      b = 2.*l*sin(alpha);    // base length of isosceles triangle
+//      d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.);   // length for InBondvector
+//      f = b/sqrt(3.);   // length for Orthvector1
+//      g = b/2.;         // length for Orthvector2
+////      LOG(3, "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", ");
+////      LOG(3, "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", "  << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g));
+//      factors[0] = d;
+//      factors[1] = f;
+//      factors[2] = 0.;
+//      FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+//      factors[1] = -0.5*f;
+//      factors[2] = g;
+//      SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+//      factors[2] = -g;
+//      ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+//
+//      // rescale each to correct BondDistance
+////      FirstOtherAtom->x.Scale(&BondRescale);
+////      SecondOtherAtom->x.Scale(&BondRescale);
+////      ThirdOtherAtom->x.Scale(&BondRescale);
+//
+//      // and relative to *origin atom
+//      *FirstOtherAtom += TopOrigin->getPosition();
+//      *SecondOtherAtom += TopOrigin->getPosition();
+//      *ThirdOtherAtom += TopOrigin->getPosition();
+//
+//      // ... and add to molecule
+//      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
+//      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
+//      AllWentWell = AllWentWell && AddAtom(ThirdOtherAtom);
+////      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
+////      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
+////      LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
+//      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
+//      Binder->Cyclic = false;
+//      Binder->Type = GraphEdge::TreeEdge;
+//      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
+//      Binder->Cyclic = false;
+//      Binder->Type = GraphEdge::TreeEdge;
+//      Binder = AddBond(BottomOrigin, ThirdOtherAtom, 1);
+//      Binder->Cyclic = false;
+//      Binder->Type = GraphEdge::TreeEdge;
+//      break;
+//    default:
+//      ELOG(1, "BondDegree does not state single, double or triple bond!");
+//      AllWentWell = false;
+//      break;
+//  }
+//
+//  return AllWentWell;
+//};
 /** Creates a copy of this molecule.

src/molecule.hpp

r2a03b0	ree4f2d
261	261	/// Add/remove atoms to/from molecule.
262	262	atom * AddCopyAtom(atom *pointer);
263		bool AddHydrogenReplacementAtom(bond::ptr Bond, atom BottomOrigin, atom TopOrigin, atom *TopReplacement, bool IsAngstroem);
	263	// bool AddHydrogenReplacementAtom(bond::ptr Bond, atom BottomOrigin, atom TopOrigin, atom *TopReplacement, bool IsAngstroem);
264	264	bond::ptr AddBond(atom first, atom second, int degree = 1);
265	265	bool hasBondStructure() const;

tests/Calculations/testsuite-calculations-1_2-dimethoxyethane.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-1_2-dimethylbenzene.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-2-methylcyclohexanone.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-N_N-dimethylacetamide.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-anthracene.at

-              r2a03b0
+              ree4f2d
         --select-molecule-by-id 0 \
         --select-molecules-atoms \
+        --correct-bonddegree \
         --fragment-molecule $FILENAME \
         --order 3 \
 …
          --fragment-resultfile ${FILENAME}_results.dat],
 , [stdout], [stderr])
 AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat | awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs($2 - energy) > 1e-5) exit(1)}'], 0)
+AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat | awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
 AT_CLEANUP

tests/Calculations/testsuite-calculations-benzene.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-cholesterol.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-coronene.at

-              r2a03b0
+              ree4f2d
         --select-molecule-by-id 0 \
         --select-molecules-atoms \
+        --correct-bonddegree \
         --fragment-molecule $FILENAME \
         --order 3 \
 …
          --fragment-resultfile ${FILENAME}_results.dat],
 , [stdout], [stderr])
 AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat | awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs($2 - energy) > 1e-5) exit(1)}'], 0)
+AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat | awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
 AT_CLEANUP

tests/Calculations/testsuite-calculations-cycloheptane.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-dimethyl_bromomalonate.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-glucose.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-heptan.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-isoleucine.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-naphthalene.at

-              r2a03b0
+              ree4f2d
         --select-molecule-by-id 0 \
         --select-molecules-atoms \
+        --correct-bonddegree \
         --fragment-molecule $FILENAME \
         --order 3 \
 …
          --fragment-resultfile ${FILENAME}_results.dat],
 , [stdout], [stderr])
 AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat | awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs($2 - energy) > 1e-5) exit(1)}'], 0)
+AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat | awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
 AT_CLEANUP

tests/Calculations/testsuite-calculations-neohexane.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-phenanthrene.at

-              r2a03b0
+              ree4f2d
         --select-molecule-by-id 0 \
         --select-molecules-atoms \
+        --correct-bonddegree \
         --fragment-molecule $FILENAME \
         --order 3 \
 …
          --fragment-resultfile ${FILENAME}_results.dat],
 , [stdout], [stderr])
 AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat | awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs($2 - energy) > 1e-5) exit(1)}'], 0)
+AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat | awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
 AT_CLEANUP

tests/Calculations/testsuite-calculations-proline.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-putrescine.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Calculations/testsuite-calculations-tartaric_acid.at

r2a03b0	ree4f2d
45	45	--fragment-resultfile ${FILENAME}_results.dat],
46	46	0, [stdout], [stderr])
47		AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(~~$2 -~~ energy) > 1e-5) exit(1)}'], 0)
	47	AT_CHECK([tail -n 1 ${FILENAME}_Energy.dat \| awk -v energy=$ENERGY 'function abs(x){return ((x < 0.0) ? -x : x)} {if (abs(($2 - energy)/energy) > 1e-5) exit(1)}'], 0)
48	48
49	49	AT_CLEANUP

tests/Fragmentations/testsuite.at

r2a03b0	ree4f2d
42	42
43	43	# Use colored output with new-enough Autotest.
44		#m4_ifdef([AT_COLOR_TESTS], [AT_COLOR_TESTS])
	44	m4_ifdef([AT_COLOR_TESTS], [AT_COLOR_TESTS])
45	45
46	46	# fragmentation of 1_2-dimethoxyethane

tests/Makefile.am

-              r2a03b0
+              ree4f2d
 extracheck:
         cd Calculations && make extracheck
+        cd Calculations && $(MAKE) extracheck
 installextracheck:
         cd Calculations && make installextracheck
+        cd Calculations && $(MAKE) installextracheck

Context Navigation

Changes in / [2a03b0:ee4f2d]

Legend:

Download in other formats: