Changeset c39675

Timestamp:

Apr 8, 2013, 11:56:08 AM (12 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

7cdf58

Parents:

6cabaac

git-author:

Frederik Heber <heber@…> (03/03/13 21:15:14)

git-committer:

Frederik Heber <heber@…> (04/08/13 11:56:08)

Message:

Added saturate() to SaturatedFragment.

• basically, this is AddHydrogenReplacementAtom() but stripped down a bit as we use hydrogens from the pool and as atoms for the fragment need not be copied anymore, we don't need all the father mambo jumbo.
• SaturatedFragment additionally now needs to know about how hydrogens are treated and whether we actually saturate.
• SaturatedFragmentUnitTest requires libLinearAlgebra.

Location:

src/Fragmentation/Exporters

Files:

• ExportGraph.cpp (modified) (1 diff)
• SaturatedFragment.cpp (modified) (3 diffs)
• SaturatedFragment.hpp (modified) (4 diffs)
• unittests/Makefile.am (modified) (2 diffs)
• unittests/SaturatedFragmentUnitTest.cpp (modified) (2 diffs)

src/Fragmentation/Exporters/ExportGraph.cpp ¶

@@ -114 +114 @@
 {
   // create the saturation which adds itself to KeySetsInUse
-  SaturatedFragment_ptr _ptr(new SaturatedFragment(_set, KeySetsInUse, hydrogens));
+  SaturatedFragment_ptr _ptr(
+      new SaturatedFragment(
+          _set,
+          KeySetsInUse,
+          hydrogens,
+          treatment,
+          saturation)
+  );
   // and return
   return _ptr;

src/Fragmentation/Exporters/SaturatedFragment.cpp ¶

@@ -37 +37 @@
 #include "SaturatedFragment.hpp"
 
+#include <cmath>
+#include <iostream>
+
 #include "CodePatterns/Assert.hpp"
-
+#include "CodePatterns/Log.hpp"
+
+#include "LinearAlgebra/Exceptions.hpp"
+#include "LinearAlgebra/Plane.hpp"
+#include "LinearAlgebra/RealSpaceMatrix.hpp"
+#include "LinearAlgebra/Vector.hpp"
+
+#include "Atom/atom.hpp"
+#include "Bond/bond.hpp"
+#include "config.hpp"
+#include "Descriptors/AtomIdDescriptor.hpp"
 #include "Fragmentation/Exporters/HydrogenPool.hpp"
+#include "Fragmentation/HydrogenSaturation_enum.hpp"
+#include "Graph/BondGraph.hpp"
+#include "World.hpp"
 
 SaturatedFragment::SaturatedFragment(
     const KeySet &_set,
     KeySetsInUse_t &_container,
-    HydrogenPool &_hydrogens) :
+    HydrogenPool &_hydrogens,
+    const enum HydrogenTreatment _treatment,
+    const enum HydrogenSaturation _saturation) :
   container(_container),
   set(_set),
   hydrogens(_hydrogens),
-  FullMolecule(set)
+  FullMolecule(set),
+  treatment(_treatment),
+  saturation(_saturation)
 {
   // add to in-use container
@@ -56 +76 @@
   container.insert(set);
 
-  // so far, we don't saturate anything
-  SaturationHydrogens.clear();
+  // prepare saturation hydrogens
+  saturate();
 }
 
@@ -77 +97 @@
   container.erase(iter);
 }
+
+void SaturatedFragment::saturate()
+{
+  // gather all atoms in a vector
+  std::vector<atom *> 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 "
+        +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();
+      ++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 in set
+      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)->BondDegree);
+////            LOG(3, output.str());
+//          //NumBonds[(*iter)->getNr()]++;
+//        } else {
+////            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
+//        }
+//        LonelyFlag = false;
+      } 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);
+        }
+//        } 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)->BondDegree);
+//        }
+        //NumBonds[(*iter)->getNr()] += Binder->BondDegree;
+      }
+    }
+  }
+}
+
+bool SaturatedFragment::OutputConfig(
+    std::ostream &out,
+    const ParserTypes _type) const
+{
+  // gather all atoms in a vector
+  std::vector<atom *> 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 "
+        +toString(*iter)+" is unknown to World.");
+    atoms.push_back(Walker);
+  }
+
+  // TODO: ScanForPeriodicCorrection() is missing so far!
+  // note however that this is not straight-forward for the following reasons:
+  // - we do not copy all atoms anymore, hence we are forced to shift the real
+  //   atoms hither and back again
+  // - we use a long-range potential that supports periodic boundary conditions.
+  //   Hence, there we would like the original configuration (split through the
+  //   the periodic boundaries). Otherwise, we would have to shift (and probably
+  //   interpolate) the potential with OBCs applying.
+
+  // list atoms in fragment for debugging
+  {
+    std::stringstream output;
+    output << "INFO: Contained atoms: ";
+    for (std::vector<atom *>::const_iterator iter = atoms.begin();
+        iter != atoms.end(); ++iter) {
+      output << (*iter)->getName() << " ";
+    }
+    LOG(3, output.str());
+  }
+
+  // store to stream via FragmentParser
+  const bool intermediateResult =
+      FormatParserStorage::getInstance().save(
+          out,
+          FormatParserStorage::getInstance().getSuffixFromType(_type),
+          atoms);
+
+  return intermediateResult;
+}
+
+atom * const SaturatedFragment::getHydrogenReplacement(atom * const replacement)
+{
+  atom * const _atom = hydrogens.leaseHydrogen();    // new atom
+  _atom->setAtomicVelocity(replacement->getAtomicVelocity()); // copy velocity
+  _atom->setFixedIon(replacement->getFixedIon());
+  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
+  _atom->father = replacement;
+  SaturationHydrogens.insert(_atom->getId());
+  return _atom;
+}
+
+bool SaturatedFragment::AddHydrogenReplacementAtom(
+    bond::ptr TopBond,
+    atom *Origin,
+    atom *Replacement,
+    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 = Replacement->getPosition() - Origin->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 Replacement or Origin being on the wrong side!
+  const BondGraph * const BG = World::getInstance().getBondGraph();
+  const range<double> MinMaxBondDistance(
+      BG->getMinMaxDistance(Origin,Replacement));
+  if (!MinMaxBondDistance.isInRange(bondlength)) {
+//    LOG(4, "InBondvector is: " << InBondvector << ".");
+    Orthovector1.Zero();
+    for (int i=NDIM;i--;) {
+      l = Replacement->at(i) - Origin->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(Origin->getType() != NULL,
+      "SaturatedFragment::AddHydrogenReplacementAtom() - element of Origin is not given.");
+  BondRescale = Origin->getType()->getHBondDistance(TopBond->BondDegree-1);
+  if (BondRescale == -1) {
+    ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->BondDegree << "!");
+    return false;
+    BondRescale = bondlength;
+  } else {
+    if (!IsAngstroem)
+      BondRescale /= (1.*AtomicLengthToAngstroem);
+  }
+
+  // discern single, double and triple bonds
+  switch(TopBond->BondDegree) {
+    case 1:
+      // check whether replacement has been an excluded hydrogen
+      if (Replacement->getType()->getAtomicNumber() == HydrogenPool::HYDROGEN) { // neither rescale nor replace if it's already hydrogen
+        FirstOtherAtom = Replacement;
+        BondRescale = bondlength;
+      } else {
+        FirstOtherAtom = getHydrogenReplacement(Replacement);
+        InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
+        FirstOtherAtom->setPosition(Origin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
+      }
+      FullMolecule.insert(FirstOtherAtom->getId());
+//      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
+      break;
+    case 2:
+      {
+        // determine two other bonds (warning if there are more than two other) plus valence sanity check
+        const BondList& ListOfBonds = Origin->getListOfBonds();
+        for (BondList::const_iterator Runner = ListOfBonds.begin();
+            Runner != ListOfBonds.end();
+            ++Runner) {
+          if ((*Runner) != TopBond) {
+            if (FirstBond == NULL) {
+              FirstBond = (*Runner);
+              FirstOtherAtom = (*Runner)->GetOtherAtom(Origin);
+            } else if (SecondBond == NULL) {
+              SecondBond = (*Runner);
+              SecondOtherAtom = (*Runner)->GetOtherAtom(Origin);
+            } else {
+              ELOG(2, "Detected more than four bonds for atom " << Origin->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 = Replacement;
+      }
+      if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
+//        LOG(3, "Regarding the double bond (" << Origin->Name << "<->" << Replacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << Origin->Name << " to determine orthogonal plane.");
+
+        // determine the plane of these two with the *origin
+        try {
+          Orthovector1 = Plane(Origin->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 = getHydrogenReplacement(Replacement);
+      SecondOtherAtom = getHydrogenReplacement(Replacement);
+      FullMolecule.insert(FirstOtherAtom->getId());
+      FullMolecule.insert(SecondOtherAtom->getId());
+      bondangle = Origin->getType()->getHBondAngle(1);
+      if (bondangle == -1) {
+        ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->BondDegree << "!");
+        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 += Origin->getPosition();
+      *SecondOtherAtom += Origin->getPosition();
+      // ... and add to molecule
+//      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
+//      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
+      break;
+    case 3:
+      // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
+      FirstOtherAtom = getHydrogenReplacement(Replacement);
+      SecondOtherAtom = getHydrogenReplacement(Replacement);
+      ThirdOtherAtom = getHydrogenReplacement(Replacement);
+      FullMolecule.insert(FirstOtherAtom->getId());
+      FullMolecule.insert(SecondOtherAtom->getId());
+      FullMolecule.insert(ThirdOtherAtom->getId());
+
+      // 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 = (Origin->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 += Origin->getPosition();
+      *SecondOtherAtom += Origin->getPosition();
+      *ThirdOtherAtom += Origin->getPosition();
+
+      // ... and add to molecule
+//      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
+//      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
+//      LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
+      break;
+    default:
+      ELOG(1, "BondDegree does not state single, double or triple bond!");
+      AllWentWell = false;
+      break;
+  }
+
+  return AllWentWell;
+};

src/Fragmentation/Exporters/SaturatedFragment.hpp ¶

@@ -14 +14 @@
 #endif
 
+#include <iosfwd>
 #include <set>
+#include <string>
 
+#include "Bond/bond.hpp"
 #include "Fragmentation/KeySet.hpp"
+#include "Fragmentation/HydrogenSaturation_enum.hpp"
+#include "Parser/FormatParserStorage.hpp"
 
+class atom;
 class HydrogenPool;
 
@@ -46 +52 @@
       const KeySet &_set,
       KeySetsInUse_t &_container,
-      HydrogenPool &_hydrogens);
+      HydrogenPool &_hydrogens,
+      const enum HydrogenTreatment _treatment,
+      const enum HydrogenSaturation saturation);
 
   /** Destructor of class SaturatedFragment.
@@ -62 +70 @@
   }
 
+  /** Getter for the FullMolecule this instance is associated with.
+   *
+   * \return const ref to FullMolecule
+   */
+  const KeySet & getFullMolecule() const
+  {
+    return FullMolecule;
+  }
+
+  /** Getter for the SaturationHydrogens this instance is associated with.
+   *
+   * \return const ref to SaturationHydrogens
+   */
+  const KeySet & getSaturationHydrogens() const
+  {
+    return SaturationHydrogens;
+  }
+
+  /** Prints the config of the fragment of \a _type to \a out.
+   *
+   * \param out output stream to write to
+   * \param _type parser type to write config
+   */
+  bool OutputConfig(
+      std::ostream &out,
+      const ParserTypes _type) const;
+
+private:
+  /** Helper function to lease and bring in place saturation hydrogens.
+   *
+   */
+  void saturate();
+
+  /** Helper function to get a hydrogen replacement for a given \a replacement.
+   *
+   * \param replacement atom to "replace" with hydrogen in a fragment.
+   * \return ref to leased hydrogen atom
+   */
+  atom * const getHydrogenReplacement(atom * const replacement);
+
+  /** Leases and adds a Hydrogen atom in replacement for the given atom \a *partner in bond with a *origin.
+   * Here, we have to distinguish between single, double or triple bonds as stated by \a BondDegree, that each demand
+   * a different scheme when adding \a *replacement atom for the given one.
+   * -# Single Bond: Simply add new atom with bond distance rescaled to typical hydrogen one
+   * -# Double Bond: Here, we need the **BondList of the \a *origin atom, by scanning for the other bonds instead of
+   *    *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalvector().
+   *    The orthonormal vector to this plane along with the vector in *Bond direction determines the plane the two
+   *    replacing hydrogens shall lie in. Now, all remains to do is take the usual hydrogen double bond angle for the
+   *    element of *origin and form the sin/cos admixture of both plane vectors for the new coordinates of the two
+   *    hydrogens forming this angle with *origin.
+   * -# Triple Bond: The idea is to set up a tetraoid (C1-H1-H2-H3) (however the lengths \f$b\f$ of the sides of the base
+   *    triangle formed by the to be added hydrogens are not equal to the typical bond distance \f$l\f$ but have to be
+   *    determined from the typical angle \f$\alpha\f$ for a hydrogen triple connected to the element of *origin):
+   *    We have the height \f$d\f$ as the vector in *Bond direction (from triangle C1-H1-H2).
+   *    \f[ h = l \cdot \cos{\left (\frac{\alpha}{2} \right )} \qquad b = 2l \cdot \sin{\left (\frac{\alpha}{2} \right)} \quad \rightarrow \quad d = l \cdot \sqrt{\cos^2{\left (\frac{\alpha}{2} \right)}-\frac{1}{3}\cdot\sin^2{\left (\frac{\alpha}{2}\right )}}
+   *    \f]
+   *    vector::GetNormalvector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalvector creates
+   *    the third one from the former two vectors. The latter ones form the plane of the base triangle mentioned above.
+   *    The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
+   *    the median lines in an isosceles triangle meet in the center point with a ratio 2:1.
+   *    \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2}
+   *    \f]
+   *    as the coordination of all three atoms in the coordinate system of these three vectors:
+   *    \f$\pmatrix{d & f & 0}\f$, \f$\pmatrix{d & -0.5 \cdot f & g}\f$ and \f$\pmatrix{d & -0.5 \cdot f & -g}\f$.
+   *
+   * \param TopBond pointer to bond between \a *origin and \a *replacement
+   * \param Origin atom that is actually contained in the fragment
+   * \param Replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
+   * \param isAngstroem whether the coordination of the given atoms is in AtomicLength (false) or Angstrom(true)
+   * \return number of atoms added, if < bond::BondDegree then something went wrong
+   */
+  bool AddHydrogenReplacementAtom(
+      bond::ptr TopBond,
+      atom *Origin,
+      atom *Replacement,
+      bool IsAngstroem);
+
 private:
   //!> container to mark ourselves RAII-style
@@ -73 +158 @@
   //!> key set containing the ids of all hydrogens added for saturation
   KeySet SaturationHydrogens;
+  //!> whether hydrogens are generally contained in set or not
+  const enum HydrogenTreatment treatment;
+  //!> whether to actually saturate or not
+  const enum HydrogenSaturation saturation;
 };
 

src/Fragmentation/Exporters/unittests/Makefile.am ¶

@@ -27 +27 @@
         ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.hpp
 HydrogenPoolUnitTest_LDADD = \
-        ${FRAGMENTATIONEXPORTERSLIBS} \
-        ../libMolecuilderUI.la
+        ../libMolecuilderUI.la \
+        ${FRAGMENTATIONEXPORTERSLIBS}
         
 SaturatedFragmentUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
@@ -34 +34 @@
         ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.hpp
 SaturatedFragmentUnitTest_LDADD = \
-        ${FRAGMENTATIONEXPORTERSLIBS} \
-        ../libMolecuilderUI.la
+        ../libMolecuilderUI.la  \
+        $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
+        ${FRAGMENTATIONEXPORTERSLIBS}
         
 

src/Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp ¶

@@ -49 +49 @@
 #include "CodePatterns/Assert.hpp"
 
+#include "Fragmentation/HydrogenSaturation_enum.hpp"
 
 #ifdef HAVE_TESTRUNNER
@@ -68 +69 @@
 
   set = new KeySet;
-  fragment = new SaturatedFragment(*set, KeySetsInUse, hydrogens);
+  fragment = new SaturatedFragment(*set, KeySetsInUse, hydrogens, ExcludeHydrogen, DoSaturate);
 }