Changes in / [ee4f2d:2a03b0]


Ignore:
Files:
2 added
11 deleted
49 edited

Legend:

Unmodified
Added
Removed

  • LinearAlgebra/src/LinearAlgebra/Makefile.am

    ree4f2d r2a03b0  
    44INCLUDES = -I$(top_srcdir)/src
    55
    6 AM_LDFLAGS = -ldl
     6AM_LDFLAGS = ${CodePatterns_LIBS} -ldl
    77AM_CPPFLAGS = $(BOOST_CPPFLAGS) ${CodePatterns_CFLAGS}
    88
    … …  
    5252libLinearAlgebra_la_LIBADD = \
    5353        $(BOOST_SERIALIZATION_LDFLAGS) $(BOOST_SERIALIZATION_LIBS) \
    54         ${CodePatterns_LIBS} \
    5554        $(GSL_LIBS)
    5655nobase_libLinearAlgebra_la_include_HEADERS = ${LINALGHEADER}

  • LinearAlgebra/src/LinearAlgebra/Plane.cpp

    ree4f2d r2a03b0  
    5757  Vector x1 = y1 - y2;
    5858  Vector x2 = y3 - y2;
    59   if ((x1.NormSquared() <= LINALG_MYEPSILON())) {
     59  if ((x1.Norm() <= LINALG_MYEPSILON())) {
    6060    throw LinearDependenceException() << LinearAlgebraVectorPair( make_pair(&y1, &y2) );
    6161  }
    62   if ((x2.NormSquared() <= LINALG_MYEPSILON())) {
     62  if ((x2.Norm() <= LINALG_MYEPSILON())) {
    6363    throw LinearDependenceException() << LinearAlgebraVectorPair( make_pair(&y2, &y3) );
    6464  }
    65   const Vector lineardependent = x1.Projection(x2) - x1;
    66   if((lineardependent.NormSquared() <= LINALG_MYEPSILON())) {
     65  if((fabs(x1.Angle(x2)) <= LINALG_MYEPSILON())) {
    6766    throw LinearDependenceException() << LinearAlgebraVectorPair( make_pair(&x1, &x2) );
    6867  }

  • LinearAlgebra/src/LinearAlgebra/Vector.cpp

    ree4f2d r2a03b0  
    290290
    291291/** Checks whether vector has all components zero.
    292  * @param epsilon numerical tolerance for equality
    293292 * @return true - vector is zero, false - vector is not
    294293 */
    295 bool Vector::IsZero(const double _epsilon) const
    296 {
    297   return (fabs(at(0))+fabs(at(1))+fabs(at(2)) < _epsilon);
     294bool Vector::IsZero() const
     295{
     296  return (fabs(at(0))+fabs(at(1))+fabs(at(2)) < LINALG_MYEPSILON());
    298297};
    299298
    300299/** Checks whether vector has length of 1.
    301  * @param epsilon numerical tolerance for equality
    302300 * @return true - vector is normalized, false - vector is not
    303301 */
    304 bool Vector::IsOne(const double _epsilon) const
    305 {
    306   return (fabs(Norm() - 1.) < _epsilon);
     302bool Vector::IsOne() const
     303{
     304  return (fabs(Norm() - 1.) < LINALG_MYEPSILON());
    307305};
    308306
    309307/** Checks whether vector is normal to \a *normal.
    310  * @param normal other supposedly normal vector
    311  * @param epsilon numerical tolerance for equality
    312308 * @return true - vector is normalized, false - vector is not
    313309 */
    314 bool Vector::IsNormalTo(const Vector &normal, const double _epsilon) const
    315 {
    316   if (ScalarProduct(normal) < _epsilon)
     310bool Vector::IsNormalTo(const Vector &normal) const
     311{
     312  if (ScalarProduct(normal) < LINALG_MYEPSILON())
    317313    return true;
    318314  else
    … …  
    320316};
    321317
    322 /** Checks whether vector is parallel to \a *parallel.
    323  * @param parallel other supposedly parallel vector
    324  * @param epsilon numerical tolerance for equality
    325  * @return true - vector is parallel/linear dependent, false - vector is not
    326  */
    327 bool Vector::IsParallelTo(const Vector &parallel, const double _epsilon) const
    328 {
    329   if (1.-fabs(ScalarProduct(parallel)/parallel.Norm()/this->Norm()) < _epsilon)
    330     return true;
    331   else
    332     return false;
    333 };
    334 
    335318/** Checks whether vector is normal to \a *normal.
    336  * @param a other vector
    337  * @param epsilon numerical tolerance for equality
    338319 * @return true - vector is normalized, false - vector is not
    339320 */
    340 bool Vector::IsEqualTo(const Vector &a, const double _epsilon) const
     321bool Vector::IsEqualTo(const Vector &a) const
    341322{
    342323  bool status = true;
    343324  for (int i=0;i<NDIM;i++) {
    344     if (fabs(at(i) - a[i]) > _epsilon)
     325    if (fabs(at(i) - a[i]) > LINALG_MYEPSILON())
    345326      status = false;
    346327  }

  • LinearAlgebra/src/LinearAlgebra/Vector.hpp

    ree4f2d r2a03b0  
    4646  double ScalarProduct(const Vector &y) const;
    4747  double Angle(const Vector &y) const;
    48   bool IsZero(const double epsilon = LINALG_MYEPSILON()) const;
    49   bool IsOne(const double epsilon = LINALG_MYEPSILON()) const;
    50   bool IsNormalTo(const Vector &normal, const double epsilon = LINALG_MYEPSILON()) const;
    51   bool IsParallelTo(const Vector &normal, const double epsilon = LINALG_MYEPSILON()) const;
    52   bool IsEqualTo(const Vector &a, const double epsilon = LINALG_MYEPSILON()) const;
     48  bool IsZero() const;
     49  bool IsOne() const;
     50  bool IsNormalTo(const Vector &normal) const;
     51  bool IsEqualTo(const Vector &a) const;
    5352
    5453  void AddVector(const Vector &y);

  • LinearAlgebra/src/LinearAlgebra/VectorContent.cpp

    ree4f2d r2a03b0  
    504504{
    505505  double factor = Norm();
    506   if (fabs(factor) > LINALG_MYEPSILON())
    507     (*this) *= 1./factor;
     506  (*this) *= 1/factor;
    508507};
    509508

  • Makefile.am

    ree4f2d r2a03b0  
    2121doc:
    2222        mkdir -p ${DX_DOCDIR}
    23         cd doc && $(MAKE) doxygen-doc
     23        cd doc && make doxygen-doc
    2424
    2525extracheck:
    26         cd tests && $(MAKE) extracheck
     26        cd tests && make extracheck
    2727installextracheck:
    28         cd tests && $(MAKE) installextracheck
     28        cd tests && make installextracheck
    2929
    3030unity:
    31         cd src && $(MAKE) unity
     31        cd src && make unity

  • src/Actions/FragmentationAction/FragmentationAction.cpp

    ree4f2d r2a03b0  
    4242#include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp"
    4343#include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"
    44 #include "Fragmentation/Exporters/SaturatedBond.hpp"
    45 #include "Fragmentation/Exporters/SaturatedFragment.hpp"
    46 #include "Fragmentation/Exporters/SaturationDistanceMaximizer.hpp"
    4744#include "Fragmentation/Fragmentation.hpp"
    4845#include "Fragmentation/Graph.hpp"
    … …  
    264261  }
    265262
    266   // create global saturation positions map
    267   SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
    268   {
    269     // go through each atom
    270     for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
    271         iter != world.endAtomSelection(); ++iter) {
    272       const atom * const _atom = iter->second;
    273 
    274       // skip hydrogens if treated special
    275       const enum HydrogenTreatment treatment =  params.HowtoTreatHydrogen.get() ? ExcludeHydrogen : IncludeHydrogen;
    276       if ((treatment == ExcludeHydrogen) && (_atom->getType()->getAtomicNumber() == 1)) {
    277         LOG(4, "DEBUG: Skipping hydrogen atom " << *_atom);
    278         continue;
    279       }
    280 
    281       // get the valence
    282       unsigned int NumberOfPoints = _atom->getElement().getNoValenceOrbitals();
    283       LOG(3, "DEBUG: There are " << NumberOfPoints
    284           << " places to fill in in total for this atom " << *_atom << ".");
    285 
    286       // check whether there are any bonds with degree larger than 1
    287       unsigned int SumOfDegrees = 0;
    288       bool PresentHigherBonds = false;
    289       const BondList &bondlist = _atom->getListOfBonds();
    290       for (BondList::const_iterator bonditer = bondlist.begin();
    291           bonditer != bondlist.end(); ++bonditer) {
    292         SumOfDegrees += (*bonditer)->getDegree();
    293         PresentHigherBonds |= (*bonditer)->getDegree() > 1;
    294       }
    295 
    296       // check whether there are alphas to maximize the hydrogens distances
    297       SaturationDistanceMaximizer::position_bins_t position_bins;
    298       {
    299         // gather all bonds and convert to SaturatedBonds
    300         SaturationDistanceMaximizer::PositionContainers_t CutBonds;
    301         for (BondList::const_iterator bonditer = bondlist.begin();
    302             bonditer != bondlist.end(); ++bonditer) {
    303           CutBonds.push_back(
    304               SaturatedBond::ptr(new SaturatedBond(*(bonditer->get()), *_atom) )
    305             );
    306         }
    307         SaturationDistanceMaximizer maximizer(CutBonds);
    308         if (PresentHigherBonds) {
    309           // then find best alphas
    310           maximizer();
    311         } else {
    312           // if no higher order bonds, we simply gather the scaled positions
    313         }
    314         position_bins = maximizer.getAllPositionBins();
    315         LOG(4, "DEBUG: Positions for atom " << *_atom << " are " << position_bins);
    316       }
    317 
    318       // convert into the desired entry in the map
    319       SaturatedFragment::SaturationsPositionsPerNeighbor_t positions_per_neighbor;
    320       {
    321         BondList::const_iterator bonditer = bondlist.begin();
    322         SaturationDistanceMaximizer::position_bins_t::const_iterator biniter =
    323             position_bins.begin();
    324 
    325         for (;bonditer != bondlist.end(); ++bonditer, ++biniter) {
    326           const atom * const OtherAtom = (*bonditer)->GetOtherAtom(_atom);
    327           std::pair<
    328               SaturatedFragment::SaturationsPositionsPerNeighbor_t::iterator,
    329               bool
    330               > inserter;
    331           // check whether we treat hydrogen special
    332           if ((treatment == ExcludeHydrogen) && (OtherAtom->getType()->getAtomicNumber() == 1)) {
    333             // if hydrogen, forget rescaled position and use original one
    334             inserter =
    335                 positions_per_neighbor.insert(
    336                     std::make_pair(
    337                         OtherAtom->getId(),
    338                         SaturatedFragment::SaturationsPositions_t(
    339                             1, OtherAtom->getPosition() - _atom->getPosition())
    340                     )
    341                 );
    342           } else {
    343             inserter =
    344                 positions_per_neighbor.insert(
    345                     std::make_pair(
    346                         OtherAtom->getId(),
    347                         SaturatedFragment::SaturationsPositions_t(
    348                             biniter->begin(),
    349                             biniter->end())
    350                     )
    351                 );
    352           }
    353           // if already pressent, add to this present list
    354           ASSERT (inserter.second,
    355               "FragmentationAction::performCall() - other atom "
    356               +toString(*OtherAtom)+" already present?");
    357         }
    358         // bonditer follows nicely
    359         ASSERT( biniter == position_bins.end(),
    360             "FragmentationAction::performCall() - biniter is out of step, it still points at bond "
    361             +toString(*biniter)+".");
    362       }
    363       // and insert
    364       globalsaturationpositions.insert(
    365           std::make_pair( _atom->getId(),
    366               positions_per_neighbor
    367           ));
    368     }
    369   }
    370 
    371263  {
    372264    const enum HydrogenSaturation saturation =  params.DoSaturation.get() ? DoSaturate : DontSaturate;
    … …  
    374266    if (params.types.get().size() != 0) {
    375267      // store molecule's fragment to file
    376       ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
     268      ExportGraph_ToFiles exporter(TotalGraph, treatment, saturation);
    377269      exporter.setPrefix(params.prefix.get());
    378270      exporter.setOutputTypes(params.types.get());
    … …  
    380272    } else {
    381273      // store molecule's fragment in FragmentJobQueue
    382       ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation, globalsaturationpositions);
     274      ExportGraph_ToJobs exporter(TotalGraph, treatment, saturation);
    383275      exporter.setLevel(params.level.get());
    384276      exporter();

  • src/Actions/FragmentationAction/StoreSaturatedFragmentAction.cpp

    ree4f2d r2a03b0  
    3939#include "CodePatterns/Log.hpp"
    4040#include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp"
    41 #include "Fragmentation/Exporters/SaturatedFragment.hpp"
    4241#include "Fragmentation/Graph.hpp"
    4342#include "World.hpp"
    … …  
    7877  // store molecule's fragment to file
    7978  {
    80     // we use an empty map here such that saturation is done locally
    81     SaturatedFragment::GlobalSaturationPositions_t globalsaturationpositions;
    82 
    8379    const enum HydrogenSaturation saturation =  params.DoSaturation.get() ? DoSaturate : DontSaturate;
    84     ExportGraph_ToFiles exporter(TotalGraph, IncludeHydrogen, saturation, globalsaturationpositions);
     80    ExportGraph_ToFiles exporter(TotalGraph, IncludeHydrogen, saturation);
    8581    exporter.setPrefix(params.prefix.get());
    8682    exporter.setOutputTypes(params.types.get());

  • src/Actions/GlobalListOfActions.hpp

    ree4f2d r2a03b0  
    2121  (Redo) \
    2222  (GraphUpdateMolecules) \
    23   (GraphCorrectBondDegree) \
    2423  (GraphCreateAdjacency) \
    2524  (GraphDepthFirstSearch) \

  • src/Actions/Makefile.am

    ree4f2d r2a03b0  
    250250
    251251GRAPHACTIONSOURCE = \
    252   Actions/GraphAction/CorrectBondDegreeAction.cpp \
    253252  Actions/GraphAction/CreateAdjacencyAction.cpp \
    254253  Actions/GraphAction/DepthFirstSearchAction.cpp \
    … …  
    257256  Actions/GraphAction/UpdateMoleculesAction.cpp
    258257GRAPHACTIONHEADER = \
    259   Actions/GraphAction/CorrectBondDegreeAction.hpp \
    260258  Actions/GraphAction/CreateAdjacencyAction.hpp \
    261259  Actions/GraphAction/DepthFirstSearchAction.hpp \
    … …  
    264262  Actions/GraphAction/UpdateMoleculesAction.hpp
    265263GRAPHACTIONDEFS = \
    266   Actions/GraphAction/CorrectBondDegreeAction.def \
    267264  Actions/GraphAction/CreateAdjacencyAction.def \
    268265  Actions/GraphAction/DepthFirstSearchAction.def \

  • src/Atom/atom_bondedparticleinfo.hpp

    ree4f2d r2a03b0  
    2323
    2424#include <list>
    25 #include <vector>
    2625
    2726/****************************************** forward declarations *****************************/
    … …  
    2928class BondedParticle;
    3029
    31 typedef std::list<bond::ptr > BondList;
     30#define BondList list<bond::ptr >
    3231
    3332/********************************************** declarations *******************************/

  • src/Element/elements_db.cpp

    ree4f2d r2a03b0  
    275275const char *HbondangleDB =\
    276276"# atomicnumber angles for single, double and triple bond (-1 no angle)\n\
    277 1       0       -1      -1\n\
    278 5       0       131.0   109.2\n\
    279 6       0       120     109.47\n\
    280 7       0       110     106.67\n\
    281 8       0       104.5   -1\n\
    282 14      0       120     109.47\n\
    283 15      0       -1      -1\n\
    284 16      0       -1      -1\n\
    285 17      0 -1 -1\n\
    286 20      0       120     109.47\n\
    287 34      0 -1 -1\n\
    288 35      0 -1 -1\n\
     2771       180     -1      -1\n\
     2785       180     131.0   109.2\n\
     2796       180     120     109.47\n\
     2807       180     110     106.67\n\
     2818       180     104.5   -1\n\
     28214      180     120     109.47\n\
     28315      180     -1      -1\n\
     28416      180     -1      -1\n\
     28517      180 -1 -1\n\
     28620      180     120     109.47\n\
     28734      180 -1 -1\n\
     28835      180 -1 -1\n\
    289289";
    290290

  • src/Fragmentation/Exporters/ExportGraph.cpp

    ree4f2d r2a03b0  
    5858  const Graph &_graph,
    5959  const enum HydrogenTreatment _treatment,
    60   const enum HydrogenSaturation _saturation,
    61   const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
     60  const enum HydrogenSaturation _saturation) :
    6261  TotalGraph(_graph),
    6362  BondFragments(World::getPointer()),
    6463  treatment(_treatment),
    6564  saturation(_saturation),
    66   globalsaturationpositions(_globalsaturationpositions),
    6765  CurrentKeySet(TotalGraph.begin())
    6866{
    … …  
    117115          hydrogens,
    118116          treatment,
    119           saturation,
    120           globalsaturationpositions)
     117          saturation)
    121118  );
    122119  // and return
    … …  
    140137}
    141138
    142 ///** Internal helper to create from each keyset a molecule
    143 // *
    144 // */
    145 //void ExportGraph::prepareMolecule()
    146 //{
    147 //  size_t count = 0;
    148 //  for(Graph::const_iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
    149 //    KeySet test = (*runner).first;
    150 //    LOG(2, "DEBUG: Fragment No." << (*runner).second.first << " with TEFactor "
    151 //        << (*runner).second.second << ".");
    152 //    BondFragments.insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
    153 //    ++count;
    154 //  }
    155 //  LOG(1, "INFO: " << count << "/" << BondFragments.ListOfMolecules.size()
    156 //      << " fragments generated from the keysets.");
    157 //}
    158 //
    159 ///** Stores a fragment from \a KeySet into \a molecule.
    160 // * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
    161 // * molecule and adds missing hydrogen where bonds were cut.
    162 // * \param &Leaflet pointer to KeySet structure
    163 // * \param IsAngstroem whether we have Ansgtroem or bohrradius
    164 // * \return pointer to constructed molecule
    165 // */
    166 //molecule * ExportGraph::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
    167 //{
    168 //  Info info(__func__);
    169 //  ListOfLocalAtoms_t SonList;
    170 //  molecule *Leaf = World::getInstance().createMolecule();
    171 //
    172 //  StoreFragmentFromKeySet_Init(Leaf, Leaflet, SonList);
    173 //  // create the bonds between all: Make it an induced subgraph and add hydrogen
    174 ////  LOG(2, "Creating bonds from father graph (i.e. induced subgraph creation).");
    175 //  CreateInducedSubgraphOfFragment(Leaf, SonList, IsAngstroem);
    176 //
    177 //  //Leaflet->Leaf->ScanForPeriodicCorrection(out);
    178 //  return Leaf;
    179 //}
    180 //
    181 ///** Initializes some value for putting fragment of \a *mol into \a *Leaf.
    182 // * \param *Leaf fragment molecule
    183 // * \param &Leaflet pointer to KeySet structure
    184 // * \param SonList calloc'd list which atom of \a *Leaf is a son of which atom in \a *mol
    185 // * \return number of atoms in fragment
    186 // */
    187 //int ExportGraph::StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList)
    188 //{
    189 //  atom *FatherOfRunner = NULL;
    190 //
    191 //  // first create the minimal set of atoms from the KeySet
    192 //  World &world = World::getInstance();
    193 //  int size = 0;
    194 //  for(KeySet::const_iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
    195 //    FatherOfRunner = world.getAtom(AtomById(*runner));  // find the id
    196 //    SonList.insert( std::make_pair(FatherOfRunner->getNr(), Leaf->AddCopyAtom(FatherOfRunner) ) );
    197 //    size++;
    198 //  }
    199 //  return size;
    200 //}
    201 //
    202 ///** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
    203 // * \param *Leaf fragment molecule
    204 // * \param IsAngstroem whether we have Ansgtroem or bohrradius
    205 // * \param SonList list which atom of \a *Leaf is another atom's son
    206 // */
    207 //void ExportGraph::CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem)
    208 //{
    209 //  bool LonelyFlag = false;
    210 //  atom *OtherFather = NULL;
    211 //  atom *FatherOfRunner = NULL;
    212 //
    213 //  // we increment the iter just before skipping the hydrogen
    214 //  // as we use AddBond, we cannot have a const_iterator here
    215 //  for (molecule::iterator iter = Leaf->begin(); iter != Leaf->end();) {
    216 //    LonelyFlag = true;
    217 //    FatherOfRunner = (*iter)->father;
    218 //    ASSERT(FatherOfRunner,"Atom without father found");
    219 //    if (SonList.find(FatherOfRunner->getNr()) != SonList.end())  {  // check if this, our father, is present in list
    220 //      // create all bonds
    221 //      const BondList& ListOfBonds = FatherOfRunner->getListOfBonds();
    222 //      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
    223 //          BondRunner != ListOfBonds.end();
    224 //          ++BondRunner) {
    225 //        OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
    226 //        if (SonList.find(OtherFather->getNr()) != SonList.end()) {
    227 ////          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
    228 ////              << " is bound to " << *OtherFather << ", whose son is "
    229 ////              << *SonList[OtherFather->getNr()] << ".");
    230 //          if (OtherFather->getNr() > FatherOfRunner->getNr()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
    231 //            std::stringstream output;
    232 ////            output << "ACCEPT: Adding Bond: "
    233 //            output << Leaf->AddBond((*iter), SonList[OtherFather->getNr()], (*BondRunner)->getDegree());
    234 ////            LOG(3, output.str());
    235 //            //NumBonds[(*iter)->getNr()]++;
    236 //          } else {
    237 ////            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
    238 //          }
    239 //          LonelyFlag = false;
    240 //        } else {
    241 ////          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
    242 ////              << " is bound to " << *OtherFather << ", who has no son in this fragment molecule.");
    243 //          if (saturation == DoSaturate) {
    244 ////          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
    245 //            if (!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
    246 //              exit(1);
    247 //          } else if ((treatment == ExcludeHydrogen) && (OtherFather->getElementNo() == (atomicNumber_t)1)) {
    248 //            // just copy the atom if it's a hydrogen
    249 //            atom * const OtherWalker = Leaf->AddCopyAtom(OtherFather);
    250 //            Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
    251 //          }
    252 //          //NumBonds[(*iter)->getNr()] += Binder->getDegree();
    253 //        }
    254 //      }
    255 //    } else {
    256 //      ELOG(1, "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->getNr()] << "!");
    257 //    }
    258 //    if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
    259 //      LOG(0, **iter << "has got bonds only to hydrogens!");
    260 //    }
    261 //    ++iter;
    262 //    if (saturation == DoSaturate) {
    263 //      while ((iter != Leaf->end()) && ((*iter)->getType()->getAtomicNumber() == 1)){ // skip added hydrogen
    264 //        iter++;
    265 //      }
    266 //    }
    267 //  }
    268 //}
     139/** Internal helper to create from each keyset a molecule
     140 *
     141 */
     142void ExportGraph::prepareMolecule()
     143{
     144  size_t count = 0;
     145  for(Graph::const_iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
     146    KeySet test = (*runner).first;
     147    LOG(2, "DEBUG: Fragment No." << (*runner).second.first << " with TEFactor "
     148        << (*runner).second.second << ".");
     149    BondFragments.insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
     150    ++count;
     151  }
     152  LOG(1, "INFO: " << count << "/" << BondFragments.ListOfMolecules.size()
     153      << " fragments generated from the keysets.");
     154}
     155
     156/** Stores a fragment from \a KeySet into \a molecule.
     157 * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
     158 * molecule and adds missing hydrogen where bonds were cut.
     159 * \param &Leaflet pointer to KeySet structure
     160 * \param IsAngstroem whether we have Ansgtroem or bohrradius
     161 * \return pointer to constructed molecule
     162 */
     163molecule * ExportGraph::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
     164{
     165  Info info(__func__);
     166  ListOfLocalAtoms_t SonList;
     167  molecule *Leaf = World::getInstance().createMolecule();
     168
     169  StoreFragmentFromKeySet_Init(Leaf, Leaflet, SonList);
     170  // create the bonds between all: Make it an induced subgraph and add hydrogen
     171//  LOG(2, "Creating bonds from father graph (i.e. induced subgraph creation).");
     172  CreateInducedSubgraphOfFragment(Leaf, SonList, IsAngstroem);
     173
     174  //Leaflet->Leaf->ScanForPeriodicCorrection(out);
     175  return Leaf;
     176}
     177
     178/** Initializes some value for putting fragment of \a *mol into \a *Leaf.
     179 * \param *Leaf fragment molecule
     180 * \param &Leaflet pointer to KeySet structure
     181 * \param SonList calloc'd list which atom of \a *Leaf is a son of which atom in \a *mol
     182 * \return number of atoms in fragment
     183 */
     184int ExportGraph::StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList)
     185{
     186  atom *FatherOfRunner = NULL;
     187
     188  // first create the minimal set of atoms from the KeySet
     189  World &world = World::getInstance();
     190  int size = 0;
     191  for(KeySet::const_iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
     192    FatherOfRunner = world.getAtom(AtomById(*runner));  // find the id
     193    SonList.insert( std::make_pair(FatherOfRunner->getNr(), Leaf->AddCopyAtom(FatherOfRunner) ) );
     194    size++;
     195  }
     196  return size;
     197}
     198
     199/** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
     200 * \param *Leaf fragment molecule
     201 * \param IsAngstroem whether we have Ansgtroem or bohrradius
     202 * \param SonList list which atom of \a *Leaf is another atom's son
     203 */
     204void ExportGraph::CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem)
     205{
     206  bool LonelyFlag = false;
     207  atom *OtherFather = NULL;
     208  atom *FatherOfRunner = NULL;
     209
     210  // we increment the iter just before skipping the hydrogen
     211  // as we use AddBond, we cannot have a const_iterator here
     212  for (molecule::iterator iter = Leaf->begin(); iter != Leaf->end();) {
     213    LonelyFlag = true;
     214    FatherOfRunner = (*iter)->father;
     215    ASSERT(FatherOfRunner,"Atom without father found");
     216    if (SonList.find(FatherOfRunner->getNr()) != SonList.end())  {  // check if this, our father, is present in list
     217      // create all bonds
     218      const BondList& ListOfBonds = FatherOfRunner->getListOfBonds();
     219      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
     220          BondRunner != ListOfBonds.end();
     221          ++BondRunner) {
     222        OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
     223        if (SonList.find(OtherFather->getNr()) != SonList.end()) {
     224//          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
     225//              << " is bound to " << *OtherFather << ", whose son is "
     226//              << *SonList[OtherFather->getNr()] << ".");
     227          if (OtherFather->getNr() > FatherOfRunner->getNr()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
     228            std::stringstream output;
     229//            output << "ACCEPT: Adding Bond: "
     230            output << Leaf->AddBond((*iter), SonList[OtherFather->getNr()], (*BondRunner)->getDegree());
     231//            LOG(3, output.str());
     232            //NumBonds[(*iter)->getNr()]++;
     233          } else {
     234//            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
     235          }
     236          LonelyFlag = false;
     237        } else {
     238//          LOG(2, "INFO: Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->getNr()]
     239//              << " is bound to " << *OtherFather << ", who has no son in this fragment molecule.");
     240          if (saturation == DoSaturate) {
     241//          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
     242            if (!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
     243              exit(1);
     244          } else if ((treatment == ExcludeHydrogen) && (OtherFather->getElementNo() == (atomicNumber_t)1)) {
     245            // just copy the atom if it's a hydrogen
     246            atom * const OtherWalker = Leaf->AddCopyAtom(OtherFather);
     247            Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
     248          }
     249          //NumBonds[(*iter)->getNr()] += Binder->getDegree();
     250        }
     251      }
     252    } else {
     253      ELOG(1, "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->getNr()] << "!");
     254    }
     255    if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
     256      LOG(0, **iter << "has got bonds only to hydrogens!");
     257    }
     258    ++iter;
     259    if (saturation == DoSaturate) {
     260      while ((iter != Leaf->end()) && ((*iter)->getType()->getAtomicNumber() == 1)){ // skip added hydrogen
     261        iter++;
     262      }
     263    }
     264  }
     265}

  • src/Fragmentation/Exporters/ExportGraph.hpp

    ree4f2d r2a03b0  
    4141      const Graph &_graph,
    4242      const enum HydrogenTreatment _treatment,
    43       const enum HydrogenSaturation _saturation,
    44       const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
     43      const enum HydrogenSaturation _saturation);
    4544  virtual ~ExportGraph();
    4645
    … …  
    7877  void releaseFragment(SaturatedFragment_ptr &_ptr);
    7978
    80 //  void prepareMolecule();
    81 //  molecule * StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem);
    82 //  int StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList);
    83 //  void CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem);
     79  void prepareMolecule();
     80  molecule * StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem);
     81  int StoreFragmentFromKeySet_Init(molecule *Leaf, KeySet &Leaflet, ListOfLocalAtoms_t &SonList);
     82  void CreateInducedSubgraphOfFragment(molecule *Leaf, ListOfLocalAtoms_t &SonList, bool IsAngstroem);
    8483
    8584protected:
    … …  
    103102  const enum HydrogenSaturation saturation;
    104103
    105   //!> Global information over all atoms with saturation positions to be used per fragment
    106   const SaturatedFragment::GlobalSaturationPositions_t &globalsaturationpositions;
    107 
    108104private:
    109105  //!> iterator pointing at the CurrentKeySet to be exported

  • src/Fragmentation/Exporters/ExportGraph_ToFiles.cpp

    ree4f2d r2a03b0  
    5757 * @param _treatment whether to always add already present hydrogens or not
    5858 * @param _saturation whether to saturate dangling bonds with hydrogen or not
    59  * @param _globalsaturationpositions possibly empty map with global information
    60  *        where to place saturation hydrogens to fulfill consistency principle
    6159 */
    6260ExportGraph_ToFiles::ExportGraph_ToFiles(
    6361    const Graph &_graph,
    6462    const enum HydrogenTreatment _treatment,
    65     const enum HydrogenSaturation _saturation,
    66     const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
    67                 ExportGraph(_graph, _treatment, _saturation, _globalsaturationpositions)
     63    const enum HydrogenSaturation _saturation) :
     64                ExportGraph(_graph, _treatment, _saturation)
    6865{}
    6966

  • src/Fragmentation/Exporters/ExportGraph_ToFiles.hpp

    ree4f2d r2a03b0  
    3333            const Graph &_graph,
    3434            const enum HydrogenTreatment _treatment,
    35             const enum HydrogenSaturation _saturation,
    36       const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
     35            const enum HydrogenSaturation _saturation);
    3736        virtual ~ExportGraph_ToFiles();
    3837

  • src/Fragmentation/Exporters/ExportGraph_ToJobs.cpp

    ree4f2d r2a03b0  
    5959    const Graph &_graph,
    6060    const enum HydrogenTreatment _treatment,
    61     const enum HydrogenSaturation _saturation,
    62     const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
    63     ExportGraph(_graph, _treatment, _saturation,_globalsaturationpositions),
     61    const enum HydrogenSaturation _saturation) :
     62    ExportGraph(_graph, _treatment, _saturation),
    6463    level(5)
    6564{}

  • src/Fragmentation/Exporters/ExportGraph_ToJobs.hpp

    ree4f2d r2a03b0  
    3535   * \param _treatment whether hydrogen is excluded in the _graph or not
    3636   * \param _saturation whether we saturate dangling bonds or not
    37    * \param _globalsaturationpositions possibly empty map with global information
    38    *        where to place saturation hydrogens to fulfill consistency principle
    3937   */
    4038  ExportGraph_ToJobs(
    4139      const Graph &_graph,
    4240      const enum HydrogenTreatment _treatment,
    43       const enum HydrogenSaturation _saturation,
    44       const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
     41      const enum HydrogenSaturation _saturation);
    4542  virtual ~ExportGraph_ToJobs();
    4643

  • src/Fragmentation/Exporters/SaturatedFragment.cpp

    ree4f2d r2a03b0  
    6363    HydrogenPool &_hydrogens,
    6464    const enum HydrogenTreatment _treatment,
    65     const enum HydrogenSaturation _saturation,
    66     const GlobalSaturationPositions_t &_globalsaturationpositions) :
     65    const enum HydrogenSaturation _saturation) :
    6766  container(_container),
    6867  set(_set),
    … …  
    7877  container.insert(set);
    7978
    80   // prepare saturation hydrogens, either using global information
    81   // or if not given, local information (created in the function)
    82   if (_globalsaturationpositions.empty())
    83     saturate();
    84   else
    85     saturate(_globalsaturationpositions);
     79  // prepare saturation hydrogens
     80  saturate();
    8681}
    8782
    … …  
    10499}
    105100
    106 typedef std::vector<atom *> atoms_t;
    107 
    108 atoms_t gatherAllAtoms(const KeySet &_FullMolecule)
    109 {
    110   atoms_t atoms;
    111   for (KeySet::const_iterator iter = _FullMolecule.begin();
    112       iter != _FullMolecule.end();
     101void SaturatedFragment::saturate()
     102{
     103  // gather all atoms in a vector
     104  std::vector<atom *> atoms;
     105  for (KeySet::const_iterator iter = FullMolecule.begin();
     106      iter != FullMolecule.end();
    113107      ++iter) {
    114108    atom * const Walker = World::getInstance().getAtom(AtomById(*iter));
    115109    ASSERT( Walker != NULL,
    116         "gatherAllAtoms() - id "
     110        "SaturatedFragment::OutputConfig() - id "
    117111        +toString(*iter)+" is unknown to World.");
    118112    atoms.push_back(Walker);
    119113  }
    120114
    121   return atoms;
    122 }
    123 
    124 typedef std::map<atom *, BondList > CutBonds_t;
    125 
    126 CutBonds_t gatherCutBonds(
    127     const atoms_t &_atoms,
    128     const KeySet &_set,
    129     const enum HydrogenTreatment _treatment)
    130 {
    131   //  bool LonelyFlag = false;
    132   CutBonds_t CutBonds;
    133   for (atoms_t::const_iterator iter = _atoms.begin();
    134       iter != _atoms.end();
     115//  bool LonelyFlag = false;
     116  for (std::vector<atom *>::const_iterator iter = atoms.begin();
     117      iter != atoms.end();
    135118      ++iter) {
    136119    atom * const Walker = *iter;
    … …  
    142125        ++BondRunner) {
    143126      atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
    144       // if other atom is in key set or is a specially treated hydrogen
    145       if (_set.find(OtherWalker->getId()) != _set.end()) {
     127      // if in set
     128      if (set.find(OtherWalker->getId()) != set.end()) {
    146129        LOG(4, "DEBUG: Walker " << *Walker << " is bound to " << *OtherWalker << ".");
    147       } else if ((_treatment == ExcludeHydrogen)
    148           && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
    149         LOG(4, "DEBUG: Walker " << *Walker << " is bound to specially treated hydrogen " <<
    150             *OtherWalker << ".");
     130//        if (OtherWalker->getId() > Walker->getId()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
     131////          std::stringstream output;
     132////            output << "ACCEPT: Adding Bond: "
     133//          output << Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
     134////            LOG(3, output.str());
     135//          //NumBonds[(*iter)->getNr()]++;
     136//        } else {
     137////            LOG(3, "REJECY: Not adding bond, labels in wrong order.");
     138//        }
     139//        LonelyFlag = false;
    151140      } else {
    152141        LOG(4, "DEBUG: Walker " << *Walker << " is bound to "
    153142            << *OtherWalker << ", who is not in this fragment molecule.");
    154           if (CutBonds.count(Walker) == 0)
    155             CutBonds.insert( std::make_pair(Walker, BondList() ));
    156           CutBonds[Walker].push_back(*BondRunner);
     143        if (saturation == DoSaturate) {
     144//          LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
     145          if (!AddHydrogenReplacementAtom(
     146              (*BondRunner),
     147              Walker,
     148              OtherWalker,
     149              World::getInstance().getConfig()->IsAngstroem == 1))
     150            exit(1);
     151        }
     152//        } else if ((treatment == ExcludeHydrogen) && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
     153//          // just copy the atom if it's a hydrogen
     154//          atom * const OtherWalker = Leaf->AddCopyAtom(OtherWalker);
     155//          Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
     156//        }
     157        //NumBonds[(*iter)->getNr()] += Binder->getDegree();
    157158      }
    158159    }
    159160  }
    160 
    161   return CutBonds;
    162 }
    163 
    164 typedef std::vector<atomId_t> atomids_t;
    165 
    166 atomids_t gatherPresentExcludedHydrogens(
    167     const atoms_t &_atoms,
    168     const KeySet &_set,
    169     const enum HydrogenTreatment _treatment)
    170 {
    171   //  bool LonelyFlag = false;
    172   atomids_t ExcludedHydrogens;
    173   for (atoms_t::const_iterator iter = _atoms.begin();
    174       iter != _atoms.end();
    175       ++iter) {
    176     atom * const Walker = *iter;
    177 
    178     // go through all bonds
    179     const BondList& ListOfBonds = Walker->getListOfBonds();
    180     for (BondList::const_iterator BondRunner = ListOfBonds.begin();
    181         BondRunner != ListOfBonds.end();
    182         ++BondRunner) {
    183       atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
    184       // if other atom is in key set or is a specially treated hydrogen
    185       if (_set.find(OtherWalker->getId()) != _set.end()) {
    186         LOG(6, "DEBUG: OtherWalker " << *OtherWalker << " is in set already.");
    187       } else if ((_treatment == ExcludeHydrogen)
    188           && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
    189         LOG(5, "DEBUG: Adding excluded hydrogen OtherWalker " << *OtherWalker << ".");
    190         ExcludedHydrogens.push_back(OtherWalker->getId());
    191       } else {
    192         LOG(6, "DEBUG: OtherWalker " << *Walker << " is not in this fragment molecule and no hydrogen.");
    193       }
    194     }
    195   }
    196 
    197   return ExcludedHydrogens;
    198 }
    199 
    200 void SaturatedFragment::saturate()
    201 {
    202   // so far, we just have a set of keys. Hence, convert to atom refs
    203   // and gather all atoms in a vector
    204   std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
    205 
    206   // go through each atom of the fragment and gather all cut bonds in list
    207   CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
    208 
    209   // add excluded hydrogens to FullMolecule if treated specially
    210   if (treatment == ExcludeHydrogen) {
    211     atomids_t ExcludedHydrogens = gatherPresentExcludedHydrogens(atoms, set, treatment);
    212     FullMolecule.insert(ExcludedHydrogens.begin(), ExcludedHydrogens.end());
    213   }
    214 
    215   // go through all cut bonds and replace with a hydrogen
    216   for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
    217       atomiter != CutBonds.end(); ++atomiter) {
    218     atom * const Walker = atomiter->first;
    219     if (!saturateAtom(Walker, atomiter->second))
    220       exit(1);
    221   }
    222 }
    223 
    224 void SaturatedFragment::saturate(
    225     const GlobalSaturationPositions_t &_globalsaturationpositions)
    226 {
    227   // so far, we just have a set of keys. Hence, convert to atom refs
    228   // and gather all atoms in a vector
    229   std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
    230 
    231   // go through each atom of the fragment and gather all cut bonds in list
    232   CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
    233 
    234   // add excluded hydrogens to FullMolecule if treated specially
    235   if (treatment == ExcludeHydrogen) {
    236     atomids_t ExcludedHydrogens = gatherPresentExcludedHydrogens(atoms, set, treatment);
    237     FullMolecule.insert(ExcludedHydrogens.begin(), ExcludedHydrogens.end());
    238   }
    239 
    240   // go through all cut bonds and replace with a hydrogen
    241   if (saturation == DoSaturate) {
    242     for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
    243         atomiter != CutBonds.end(); ++atomiter) {
    244       atom * const Walker = atomiter->first;
    245       LOG(4, "DEBUG: We are now saturating dangling bonds of " << *Walker);
    246 
    247       // gather set of positions for this atom from global map
    248       GlobalSaturationPositions_t::const_iterator mapiter =
    249           _globalsaturationpositions.find(Walker->getId());
    250       ASSERT( mapiter != _globalsaturationpositions.end(),
    251           "SaturatedFragment::saturate() - no global information for "
    252           +toString(*Walker));
    253       const SaturationsPositionsPerNeighbor_t &saturationpositions =
    254           mapiter->second;
    255 
    256       // go through all cut bonds for this atom
    257       for (BondList::const_iterator bonditer = atomiter->second.begin();
    258           bonditer != atomiter->second.end(); ++bonditer) {
    259         atom * const OtherWalker = (*bonditer)->GetOtherAtom(Walker);
    260 
    261         // get positions from global map
    262         SaturationsPositionsPerNeighbor_t::const_iterator positionsiter =
    263             saturationpositions.find(OtherWalker->getId());
    264         ASSERT(positionsiter != saturationpositions.end(),
    265             "SaturatedFragment::saturate() - no information on bond neighbor "
    266             +toString(*OtherWalker)+" to atom "+toString(*Walker));
    267         ASSERT(!positionsiter->second.empty(),
    268             "SaturatedFragment::saturate() - no positions for saturating bond to"
    269             +toString(*OtherWalker)+" to atom "+toString(*Walker));
    270 
    271 //        // get typical bond distance from elements database
    272 //        double BondDistance = Walker->getType()->getHBondDistance(positionsiter->second.size()-1);
    273 //        if (BondDistance < 0.) {
    274 //          ELOG(2, "saturateAtoms() - no typical hydrogen bond distance of degree "
    275 //              +toString(positionsiter->second.size())+" for element "
    276 //              +toString(Walker->getType()->getName()));
    277 //          // try bond degree 1 distance
    278 //          BondDistance = Walker->getType()->getHBondDistance(1-1);
    279 //          if (BondDistance < 0.) {
    280 //            ELOG(1, "saturateAtoms() - no typical hydrogen bond distance for element "
    281 //                +toString(Walker->getType()->getName()));
    282 //            BondDistance = 1.;
    283 //          }
    284 //        }
    285 
    286         // place hydrogen at each point
    287         LOG(4, "DEBUG: Places to saturate for atom " << *OtherWalker
    288             << " are " << positionsiter->second);
    289         atom * const father = Walker;
    290         for (SaturationsPositions_t::const_iterator positer = positionsiter->second.begin();
    291             positer != positionsiter->second.end(); ++positer) {
    292           const atom& hydrogen =
    293               setHydrogenReplacement(Walker, *positer, 1., father);
    294           FullMolecule.insert(hydrogen.getId());
    295         }
    296       }
    297     }
    298   } else
    299     LOG(3, "INFO: We are not saturating cut bonds.");
    300 }
    301 
    302 const atom& SaturatedFragment::setHydrogenReplacement(
    303   const atom * const _OwnerAtom,
    304   const Vector &_position,
    305   const double _distance,
    306   atom * const _father)
    307 {
    308   atom * const _atom = hydrogens.leaseHydrogen();    // new atom
    309   _atom->setPosition( _OwnerAtom->getPosition() + _distance * _position );
    310   // always set as fixed ion (not moving during molecular dynamics simulation)
    311   _atom->setFixedIon(true);
    312   // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
    313   _atom->father = _father;
    314   SaturationHydrogens.insert(_atom->getId());
    315 
    316   return *_atom;
    317 }
    318 
    319 bool SaturatedFragment::saturateAtom(
    320     atom * const _atom,
    321     const BondList &_cutbonds)
    322 {
    323   // go through each bond and replace
    324   for (BondList::const_iterator bonditer = _cutbonds.begin();
    325       bonditer != _cutbonds.end(); ++bonditer) {
    326     atom * const OtherWalker = (*bonditer)->GetOtherAtom(_atom);
    327     if (!AddHydrogenReplacementAtom(
    328         (*bonditer),
    329         _atom,
    330         OtherWalker,
    331         World::getInstance().getConfig()->IsAngstroem == 1))
    332       return false;
    333   }
    334   return true;
    335161}
    336162
    … …  
    444270  if (BondRescale == -1) {
    445271    ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->getDegree() << "!");
    446     BondRescale = Origin->getType()->getHBondDistance(TopBond->getDegree());
    447     if (BondRescale == -1) {
    448       ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of any degree!");
    449       return false;
    450       BondRescale = bondlength;
    451     }
     272    return false;
     273    BondRescale = bondlength;
    452274  } else {
    453275    if (!IsAngstroem)

  • src/Fragmentation/Exporters/SaturatedFragment.hpp

    ree4f2d r2a03b0  
    1818#include <string>
    1919
    20 #include "Atom/atom_bondedparticleinfo.hpp"
    2120#include "Bond/bond.hpp"
    2221#include "Fragmentation/KeySet.hpp"
    2322#include "Fragmentation/HydrogenSaturation_enum.hpp"
    2423#include "Parser/FormatParserStorage.hpp"
    25 
    26 #include "LinearAlgebra/Vector.hpp"
    2724
    2825class atom;
    … …  
    4542  typedef std::set<KeySet> KeySetsInUse_t;
    4643
    47   //!> List of points giving saturation positions for a single bond neighbor
    48   typedef std::list<Vector> SaturationsPositions_t;
    49   //!> map for one atom, containing the saturation points for all its neighbors
    50   typedef std::map<int, SaturationsPositions_t> SaturationsPositionsPerNeighbor_t;
    51   //!> containing the saturation points over all desired atoms required
    52   typedef std::map<int, SaturationsPositionsPerNeighbor_t> GlobalSaturationPositions_t;
    53 
    5444  /** Constructor of SaturatedFragment requires \a set which we are tightly
    5545   * associated.
    … …  
    5848   * \param _container container to add KeySet as in-use
    5949   * \param _hydrogens pool with hydrogens for saturation
    60    * \param _globalsaturationpositions saturation positions to be used
    6150   */
    6251  SaturatedFragment(
    … …  
    6554      HydrogenPool &_hydrogens,
    6655      const enum HydrogenTreatment _treatment,
    67       const enum HydrogenSaturation saturation,
    68       const GlobalSaturationPositions_t &_globalsaturationpositions);
     56      const enum HydrogenSaturation saturation);
    6957
    7058  /** Destructor of class SaturatedFragment.
    … …  
    112100  /** Helper function to lease and bring in place saturation hydrogens.
    113101   *
    114    * Here, we use local information to calculate saturation positions.
    115    *
    116102   */
    117103  void saturate();
    118 
    119   /** Helper function to lease and bring in place saturation hydrogens.
    120    *
    121    * Here, saturation positions have to be calculated before and are fully
    122    * stored in \a _globalsaturationpositions.
    123    *
    124    * \param_globalsaturationpositions
    125    */
    126   void saturate(const GlobalSaturationPositions_t &_globalsaturationpositions);
    127 
    128   /** Replaces all cut bonds with respect to the given atom by hydrogens.
    129    *
    130    * \param _atom atom whose cut bonds to saturate
    131    * \param _cutbonds list of cut bonds for \a _atom
    132    * \return true - bonds saturated, false - something went wrong
    133    */
    134   bool saturateAtom(atom * const _atom, const BondList &_cutbonds);
    135104
    136105  /** Helper function to get a hydrogen replacement for a given \a replacement.
    … …  
    140109   */
    141110  atom * const getHydrogenReplacement(atom * const replacement);
    142 
    143   /** Sets a saturation hydrogen at the given position in place of \a _father.
    144    *
    145    * \param _OwnerAtom atom "owning" the hydrogen (i.e. the one getting saturated)
    146    * \param _position new position relative to \a _OwnerAtom
    147    * \param _distance scale factor to the distance (default 1.)
    148    * \param _father bond partner of \a _OwnerAtom that is replaced
    149    */
    150   const atom& setHydrogenReplacement(
    151       const atom * const _OwnerAtom,
    152       const Vector &_position,
    153       const double _distance,
    154       atom * const _father);
    155111
    156112  /** 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

    ree4f2d r2a03b0  
    44FRAGMENTATIONEXPORTERSSOURCES = \
    55        ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.cpp \
    6         ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp \
    7         ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.cpp
     6        ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp
    87
    98FRAGMENTATIONEXPORTERSTESTSHEADERS = \
    109        ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.hpp \
    11         ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.hpp \
    12         ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.hpp
     10        ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.hpp
    1311
    1412FRAGMENTATIONEXPORTERSTESTS = \
    1513        HydrogenPoolUnitTest \
    16         SaturatedFragmentUnitTest \
    17         SaturationDistanceMaximizerUnitTest
     14        SaturatedFragmentUnitTest
    1815
    1916TESTS += $(FRAGMENTATIONEXPORTERSTESTS)
    … …  
    4037        $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
    4138        ${FRAGMENTATIONEXPORTERSLIBS}
    42 
    43 SaturationDistanceMaximizerUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
    44         ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.cpp \
    45         ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.hpp \
    46         ../Fragmentation/Exporters/unittests/stubs/SaturatedBondStub.cpp
    47 SaturationDistanceMaximizerUnitTest_LDADD = \
    48         ../libMolecuilderUI.la  \
    49         $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
    50         ${FRAGMENTATIONEXPORTERSLIBS}
    5139       
    5240

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

    ree4f2d r2a03b0  
    6969  ASSERT_DO(Assert::Throw);
    7070
    71   SaturatedFragment::GlobalSaturationPositions_t globalpositions;
    7271  set = new KeySet;
    73   fragment = new SaturatedFragment(
    74       *set,
    75       KeySetsInUse,
    76       hydrogens,
    77       ExcludeHydrogen,
    78       DoSaturate,
    79       globalpositions);
     72  fragment = new SaturatedFragment(*set, KeySetsInUse, hydrogens, ExcludeHydrogen, DoSaturate);
    8073}
    8174

  • src/Fragmentation/Makefile.am

    ree4f2d r2a03b0  
    77        Fragmentation/Exporters/ExportGraph.cpp \
    88        Fragmentation/Exporters/HydrogenPool.cpp \
    9         Fragmentation/Exporters/SaturatedBond.cpp \
    109        Fragmentation/Exporters/SaturatedFragment.cpp \
    11         Fragmentation/Exporters/SaturationDistanceMaximizer.cpp \
    1210        Fragmentation/Homology/FragmentEdge.cpp \
    1311        Fragmentation/Homology/FragmentNode.cpp \
    … …  
    3533        Fragmentation/Exporters/ExportGraph.hpp \
    3634        Fragmentation/Exporters/HydrogenPool.hpp \
    37         Fragmentation/Exporters/SaturatedBond.hpp \
    3835        Fragmentation/Exporters/SaturatedFragment.hpp \
    39         Fragmentation/Exporters/SaturationDistanceMaximizer.hpp \
    4036        Fragmentation/Homology/FragmentEdge.hpp \
    4137        Fragmentation/Homology/FragmentNode.hpp \

  • src/Fragmentation/Summation/Converter/DataConverter.hpp

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

  • src/Graph/Makefile.am

    ree4f2d r2a03b0  
    44GRAPHSOURCE = \
    55        Graph/BondGraph.cpp \
     6        Graph/BreadthFirstSearchAdd.cpp \
    67        Graph/BuildInducedSubgraph.cpp \
    78        Graph/AdjacencyList.cpp \
    … …  
    1213GRAPHHEADER = \
    1314        Graph/BondGraph.hpp \
     15        Graph/BreadthFirstSearchAdd.hpp \
    1416        Graph/BuildInducedSubgraph.hpp \
    1517        Graph/AdjacencyList.hpp \

  • src/Makefile.am

    ree4f2d r2a03b0  
    3636include RandomNumbers/Makefile.am
    3737include Shapes/Makefile.am
    38 include Tesselation/Makefile.am
    3938include UIElements/Makefile.am
    4039
    … …  
    147146  Thermostats/Woodcock.hpp
    148147
     148TESSELATIONSOURCE = \
     149  Tesselation/ApproximateShapeArea.cpp \
     150  Tesselation/ApproximateShapeVolume.cpp \
     151  Tesselation/boundary.cpp \
     152  Tesselation/BoundaryLineSet.cpp \
     153  Tesselation/BoundaryPointSet.cpp \
     154  Tesselation/BoundaryPolygonSet.cpp \
     155  Tesselation/BoundaryTriangleSet.cpp \
     156  Tesselation/CandidateForTesselation.cpp \
     157  Tesselation/ellipsoid.cpp \
     158  Tesselation/tesselation.cpp \
     159  Tesselation/tesselationhelpers.cpp \
     160  Tesselation/triangleintersectionlist.cpp
     161 
     162TESSELATIONHEADER = \
     163  Tesselation/ApproximateShapeArea.hpp \
     164  Tesselation/ApproximateShapeVolume.hpp \
     165  Tesselation/boundary.hpp \
     166  Tesselation/BoundaryLineSet.hpp \
     167  Tesselation/BoundaryMaps.hpp \
     168  Tesselation/BoundaryPointSet.hpp \
     169  Tesselation/BoundaryPolygonSet.hpp \
     170  Tesselation/BoundaryTriangleSet.hpp \
     171  Tesselation/CandidateForTesselation.hpp \
     172  Tesselation/ellipsoid.hpp \
     173  Tesselation/tesselation.hpp \
     174  Tesselation/tesselationhelpers.hpp \
     175  Tesselation/triangleintersectionlist.hpp
     176
    149177MOLECUILDERSOURCE = \
    150178  ${BONDSOURCE} \
    … …  
    152180  ${DYNAMICSSOURCE} \
    153181  ${THERMOSTATSOURCE} \
     182  ${TESSELATIONSOURCE} \
    154183  Shapes/ShapeFactory.cpp \
    155184  AtomIdSet.cpp \
    … …  
    174203  ${DYNAMICSHEADER} \
    175204  ${THERMOSTATHEADER} \
     205  ${TESSELATIONHEADER} \
    176206  Shapes/ShapeFactory.hpp \
    177207  AtomIdSet.hpp \
    … …  
    201231        $(BOOST_THREAD_LDFLAGS)
    202232libMolecuilder_la_LIBADD = \
    203         libMolecuilderTesselation.la \
    204233        libMolecuilderShapes.la \
    205234        $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \

  • src/molecule.cpp

    ree4f2d r2a03b0  
    339339 * \todo double and triple bonds splitting (always use the tetraeder angle!)
    340340 */
    341 //bool molecule::AddHydrogenReplacementAtom(bond::ptr TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
    342 //{
    343 ////  Info info(__func__);
    344 //  bool AllWentWell = true;    // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
    345 //  double bondlength;  // bond length of the bond to be replaced/cut
    346 //  double bondangle;  // bond angle of the bond to be replaced/cut
    347 //  double BondRescale;   // rescale value for the hydrogen bond length
    348 //  bond::ptr FirstBond;
    349 //  bond::ptr SecondBond; // Other bonds in double bond case to determine "other" plane
    350 //  atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
    351 //  double b,l,d,f,g, alpha, factors[NDIM];    // hold temporary values in triple bond case for coordination determination
    352 //  Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
    353 //  Vector InBondvector;    // vector in direction of *Bond
    354 //  const RealSpaceMatrix &matrix =  World::getInstance().getDomain().getM();
    355 //  bond::ptr Binder;
    356 //
    357 //  // create vector in direction of bond
    358 //  InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
    359 //  bondlength = InBondvector.Norm();
    360 //
    361 //   // is greater than typical bond distance? Then we have to correct periodically
    362 //   // the problem is not the H being out of the box, but InBondvector have the wrong direction
    363 //   // due to TopReplacement or Origin being on the wrong side!
    364 //  const BondGraph * const BG = World::getInstance().getBondGraph();
    365 //  const range<double> MinMaxBondDistance(
    366 //      BG->getMinMaxDistance(TopOrigin,TopReplacement));
    367 //  if (!MinMaxBondDistance.isInRange(bondlength)) {
    368 ////    LOG(4, "InBondvector is: " << InBondvector << ".");
    369 //    Orthovector1.Zero();
    370 //    for (int i=NDIM;i--;) {
    371 //      l = TopReplacement->at(i) - TopOrigin->at(i);
    372 //      if (fabs(l) > MinMaxBondDistance.last) { // is component greater than bond distance (check against min not useful here)
    373 //        Orthovector1[i] = (l < 0) ? -1. : +1.;
    374 //      } // (signs are correct, was tested!)
    375 //    }
    376 //    Orthovector1 *= matrix;
    377 //    InBondvector -= Orthovector1; // subtract just the additional translation
    378 //    bondlength = InBondvector.Norm();
    379 ////    LOG(4, "INFO: Corrected InBondvector is now: " << InBondvector << ".");
    380 //  } // periodic correction finished
    381 //
    382 //  InBondvector.Normalize();
    383 //  // get typical bond length and store as scale factor for later
    384 //  ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
    385 //  BondRescale = TopOrigin->getType()->getHBondDistance(TopBond->getDegree()-1);
    386 //  if (BondRescale == -1) {
    387 //    ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
    388 //    return false;
    389 //    BondRescale = bondlength;
    390 //  } else {
    391 //    if (!IsAngstroem)
    392 //      BondRescale /= (1.*AtomicLengthToAngstroem);
    393 //  }
    394 //
    395 //  // discern single, double and triple bonds
    396 //  switch(TopBond->getDegree()) {
    397 //    case 1:
    398 //      FirstOtherAtom = World::getInstance().createAtom();    // new atom
    399 //      FirstOtherAtom->setType(1);  // element is Hydrogen
    400 //      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
    401 //      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
    402 //      if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
    403 //        FirstOtherAtom->father = TopReplacement;
    404 //        BondRescale = bondlength;
    405 //      } else {
    406 //        FirstOtherAtom->father = NULL;  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
    407 //      }
    408 //      InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
    409 //      FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
    410 //      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
    411 ////      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
    412 //      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
    413 //      Binder->Cyclic = false;
    414 //      Binder->Type = GraphEdge::TreeEdge;
    415 //      break;
    416 //    case 2:
    417 //      {
    418 //        // determine two other bonds (warning if there are more than two other) plus valence sanity check
    419 //        const BondList& ListOfBonds = TopOrigin->getListOfBonds();
    420 //        for (BondList::const_iterator Runner = ListOfBonds.begin();
    421 //            Runner != ListOfBonds.end();
    422 //            ++Runner) {
    423 //          if ((*Runner) != TopBond) {
    424 //            if (FirstBond == NULL) {
    425 //              FirstBond = (*Runner);
    426 //              FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
    427 //            } else if (SecondBond == NULL) {
    428 //              SecondBond = (*Runner);
    429 //              SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
    430 //            } else {
    431 //              ELOG(2, "Detected more than four bonds for atom " << TopOrigin->getName());
    432 //            }
    433 //          }
    434 //        }
    435 //      }
    436 //      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)
    437 //        SecondBond = TopBond;
    438 //        SecondOtherAtom = TopReplacement;
    439 //      }
    440 //      if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
    441 ////        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.");
    442 //
    443 //        // determine the plane of these two with the *origin
    444 //        try {
    445 //          Orthovector1 = Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
    446 //        }
    447 //        catch(LinearDependenceException &excp){
    448 //          LOG(0, boost::diagnostic_information(excp));
    449 //          // TODO: figure out what to do with the Orthovector in this case
    450 //          AllWentWell = false;
    451 //        }
    452 //      } else {
    453 //        Orthovector1.GetOneNormalVector(InBondvector);
    454 //      }
    455 //      //LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
    456 //      // orthogonal vector and bond vector between origin and replacement form the new plane
    457 //      Orthovector1.MakeNormalTo(InBondvector);
    458 //      Orthovector1.Normalize();
    459 //      //LOG(3, "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << ".");
    460 //
    461 //      // create the two Hydrogens ...
    462 //      FirstOtherAtom = World::getInstance().createAtom();
    463 //      SecondOtherAtom = World::getInstance().createAtom();
    464 //      FirstOtherAtom->setType(1);
    465 //      SecondOtherAtom->setType(1);
    466 //      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
    467 //      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
    468 //      SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
    469 //      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
    470 //      FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
    471 //      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
    472 //      bondangle = TopOrigin->getType()->getHBondAngle(1);
    473 //      if (bondangle == -1) {
    474 //        ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
    475 //        return false;
    476 //        bondangle = 0;
    477 //      }
    478 //      bondangle *= M_PI/180./2.;
    479 ////      LOG(3, "INFO: ReScaleCheck: InBondvector " << InBondvector << ", " << Orthovector1 << ".");
    480 ////      LOG(3, "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle));
    481 //      FirstOtherAtom->Zero();
    482 //      SecondOtherAtom->Zero();
    483 //      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
    484 //        FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
    485 //        SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
    486 //      }
    487 //      FirstOtherAtom->Scale(BondRescale);  // rescale by correct BondDistance
    488 //      SecondOtherAtom->Scale(BondRescale);
    489 //      //LOG(3, "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << ".");
    490 //      *FirstOtherAtom += TopOrigin->getPosition();
    491 //      *SecondOtherAtom += TopOrigin->getPosition();
    492 //      // ... and add to molecule
    493 //      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
    494 //      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
    495 ////      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
    496 ////      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
    497 //      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
    498 //      Binder->Cyclic = false;
    499 //      Binder->Type = GraphEdge::TreeEdge;
    500 //      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
    501 //      Binder->Cyclic = false;
    502 //      Binder->Type = GraphEdge::TreeEdge;
    503 //      break;
    504 //    case 3:
    505 //      // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
    506 //      FirstOtherAtom = World::getInstance().createAtom();
    507 //      SecondOtherAtom = World::getInstance().createAtom();
    508 //      ThirdOtherAtom = World::getInstance().createAtom();
    509 //      FirstOtherAtom->setType(1);
    510 //      SecondOtherAtom->setType(1);
    511 //      ThirdOtherAtom->setType(1);
    512 //      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
    513 //      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
    514 //      SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
    515 //      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
    516 //      ThirdOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
    517 //      ThirdOtherAtom->setFixedIon(TopReplacement->getFixedIon());
    518 //      FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
    519 //      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
    520 //      ThirdOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
    521 //
    522 //      // we need to vectors orthonormal the InBondvector
    523 //      AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
    524 ////      LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
    525 //      try{
    526 //        Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
    527 //      }
    528 //      catch(LinearDependenceException &excp) {
    529 //        LOG(0, boost::diagnostic_information(excp));
    530 //        AllWentWell = false;
    531 //      }
    532 ////      LOG(3, "INFO: Orthovector2: " << Orthovector2 << ".")
    533 //
    534 //      // create correct coordination for the three atoms
    535 //      alpha = (TopOrigin->getType()->getHBondAngle(2))/180.*M_PI/2.;  // retrieve triple bond angle from database
    536 //      l = BondRescale;        // desired bond length
    537 //      b = 2.*l*sin(alpha);    // base length of isosceles triangle
    538 //      d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.);   // length for InBondvector
    539 //      f = b/sqrt(3.);   // length for Orthvector1
    540 //      g = b/2.;         // length for Orthvector2
    541 ////      LOG(3, "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", ");
    542 ////      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));
    543 //      factors[0] = d;
    544 //      factors[1] = f;
    545 //      factors[2] = 0.;
    546 //      FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
    547 //      factors[1] = -0.5*f;
    548 //      factors[2] = g;
    549 //      SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
    550 //      factors[2] = -g;
    551 //      ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
    552 //
    553 //      // rescale each to correct BondDistance
    554 ////      FirstOtherAtom->x.Scale(&BondRescale);
    555 ////      SecondOtherAtom->x.Scale(&BondRescale);
    556 ////      ThirdOtherAtom->x.Scale(&BondRescale);
    557 //
    558 //      // and relative to *origin atom
    559 //      *FirstOtherAtom += TopOrigin->getPosition();
    560 //      *SecondOtherAtom += TopOrigin->getPosition();
    561 //      *ThirdOtherAtom += TopOrigin->getPosition();
    562 //
    563 //      // ... and add to molecule
    564 //      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
    565 //      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
    566 //      AllWentWell = AllWentWell && AddAtom(ThirdOtherAtom);
    567 ////      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
    568 ////      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
    569 ////      LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
    570 //      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
    571 //      Binder->Cyclic = false;
    572 //      Binder->Type = GraphEdge::TreeEdge;
    573 //      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
    574 //      Binder->Cyclic = false;
    575 //      Binder->Type = GraphEdge::TreeEdge;
    576 //      Binder = AddBond(BottomOrigin, ThirdOtherAtom, 1);
    577 //      Binder->Cyclic = false;
    578 //      Binder->Type = GraphEdge::TreeEdge;
    579 //      break;
    580 //    default:
    581 //      ELOG(1, "BondDegree does not state single, double or triple bond!");
    582 //      AllWentWell = false;
    583 //      break;
    584 //  }
    585 //
    586 //  return AllWentWell;
    587 //};
     341bool molecule::AddHydrogenReplacementAtom(bond::ptr TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
     342{
     343//  Info info(__func__);
     344  bool AllWentWell = true;    // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
     345  double bondlength;  // bond length of the bond to be replaced/cut
     346  double bondangle;  // bond angle of the bond to be replaced/cut
     347  double BondRescale;   // rescale value for the hydrogen bond length
     348  bond::ptr FirstBond;
     349  bond::ptr SecondBond; // Other bonds in double bond case to determine "other" plane
     350  atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
     351  double b,l,d,f,g, alpha, factors[NDIM];    // hold temporary values in triple bond case for coordination determination
     352  Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
     353  Vector InBondvector;    // vector in direction of *Bond
     354  const RealSpaceMatrix &matrix =  World::getInstance().getDomain().getM();
     355  bond::ptr Binder;
     356
     357  // create vector in direction of bond
     358  InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
     359  bondlength = InBondvector.Norm();
     360
     361   // is greater than typical bond distance? Then we have to correct periodically
     362   // the problem is not the H being out of the box, but InBondvector have the wrong direction
     363   // due to TopReplacement or Origin being on the wrong side!
     364  const BondGraph * const BG = World::getInstance().getBondGraph();
     365  const range<double> MinMaxBondDistance(
     366      BG->getMinMaxDistance(TopOrigin,TopReplacement));
     367  if (!MinMaxBondDistance.isInRange(bondlength)) {
     368//    LOG(4, "InBondvector is: " << InBondvector << ".");
     369    Orthovector1.Zero();
     370    for (int i=NDIM;i--;) {
     371      l = TopReplacement->at(i) - TopOrigin->at(i);
     372      if (fabs(l) > MinMaxBondDistance.last) { // is component greater than bond distance (check against min not useful here)
     373        Orthovector1[i] = (l < 0) ? -1. : +1.;
     374      } // (signs are correct, was tested!)
     375    }
     376    Orthovector1 *= matrix;
     377    InBondvector -= Orthovector1; // subtract just the additional translation
     378    bondlength = InBondvector.Norm();
     379//    LOG(4, "INFO: Corrected InBondvector is now: " << InBondvector << ".");
     380  } // periodic correction finished
     381
     382  InBondvector.Normalize();
     383  // get typical bond length and store as scale factor for later
     384  ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
     385  BondRescale = TopOrigin->getType()->getHBondDistance(TopBond->getDegree()-1);
     386  if (BondRescale == -1) {
     387    ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
     388    return false;
     389    BondRescale = bondlength;
     390  } else {
     391    if (!IsAngstroem)
     392      BondRescale /= (1.*AtomicLengthToAngstroem);
     393  }
     394
     395  // discern single, double and triple bonds
     396  switch(TopBond->getDegree()) {
     397    case 1:
     398      FirstOtherAtom = World::getInstance().createAtom();    // new atom
     399      FirstOtherAtom->setType(1);  // element is Hydrogen
     400      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
     401      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
     402      if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
     403        FirstOtherAtom->father = TopReplacement;
     404        BondRescale = bondlength;
     405      } else {
     406        FirstOtherAtom->father = NULL;  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
     407      }
     408      InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
     409      FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
     410      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
     411//      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
     412      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
     413      Binder->Cyclic = false;
     414      Binder->Type = GraphEdge::TreeEdge;
     415      break;
     416    case 2:
     417      {
     418        // determine two other bonds (warning if there are more than two other) plus valence sanity check
     419        const BondList& ListOfBonds = TopOrigin->getListOfBonds();
     420        for (BondList::const_iterator Runner = ListOfBonds.begin();
     421            Runner != ListOfBonds.end();
     422            ++Runner) {
     423          if ((*Runner) != TopBond) {
     424            if (FirstBond == NULL) {
     425              FirstBond = (*Runner);
     426              FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
     427            } else if (SecondBond == NULL) {
     428              SecondBond = (*Runner);
     429              SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
     430            } else {
     431              ELOG(2, "Detected more than four bonds for atom " << TopOrigin->getName());
     432            }
     433          }
     434        }
     435      }
     436      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)
     437        SecondBond = TopBond;
     438        SecondOtherAtom = TopReplacement;
     439      }
     440      if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
     441//        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.");
     442
     443        // determine the plane of these two with the *origin
     444        try {
     445          Orthovector1 = Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
     446        }
     447        catch(LinearDependenceException &excp){
     448          LOG(0, boost::diagnostic_information(excp));
     449          // TODO: figure out what to do with the Orthovector in this case
     450          AllWentWell = false;
     451        }
     452      } else {
     453        Orthovector1.GetOneNormalVector(InBondvector);
     454      }
     455      //LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
     456      // orthogonal vector and bond vector between origin and replacement form the new plane
     457      Orthovector1.MakeNormalTo(InBondvector);
     458      Orthovector1.Normalize();
     459      //LOG(3, "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << ".");
     460
     461      // create the two Hydrogens ...
     462      FirstOtherAtom = World::getInstance().createAtom();
     463      SecondOtherAtom = World::getInstance().createAtom();
     464      FirstOtherAtom->setType(1);
     465      SecondOtherAtom->setType(1);
     466      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
     467      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
     468      SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
     469      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
     470      FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
     471      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
     472      bondangle = TopOrigin->getType()->getHBondAngle(1);
     473      if (bondangle == -1) {
     474        ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
     475        return false;
     476        bondangle = 0;
     477      }
     478      bondangle *= M_PI/180./2.;
     479//      LOG(3, "INFO: ReScaleCheck: InBondvector " << InBondvector << ", " << Orthovector1 << ".");
     480//      LOG(3, "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle));
     481      FirstOtherAtom->Zero();
     482      SecondOtherAtom->Zero();
     483      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
     484        FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
     485        SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
     486      }
     487      FirstOtherAtom->Scale(BondRescale);  // rescale by correct BondDistance
     488      SecondOtherAtom->Scale(BondRescale);
     489      //LOG(3, "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << ".");
     490      *FirstOtherAtom += TopOrigin->getPosition();
     491      *SecondOtherAtom += TopOrigin->getPosition();
     492      // ... and add to molecule
     493      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
     494      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
     495//      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
     496//      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
     497      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
     498      Binder->Cyclic = false;
     499      Binder->Type = GraphEdge::TreeEdge;
     500      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
     501      Binder->Cyclic = false;
     502      Binder->Type = GraphEdge::TreeEdge;
     503      break;
     504    case 3:
     505      // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
     506      FirstOtherAtom = World::getInstance().createAtom();
     507      SecondOtherAtom = World::getInstance().createAtom();
     508      ThirdOtherAtom = World::getInstance().createAtom();
     509      FirstOtherAtom->setType(1);
     510      SecondOtherAtom->setType(1);
     511      ThirdOtherAtom->setType(1);
     512      FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
     513      FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
     514      SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
     515      SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
     516      ThirdOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
     517      ThirdOtherAtom->setFixedIon(TopReplacement->getFixedIon());
     518      FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
     519      SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
     520      ThirdOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
     521
     522      // we need to vectors orthonormal the InBondvector
     523      AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
     524//      LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
     525      try{
     526        Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
     527      }
     528      catch(LinearDependenceException &excp) {
     529        LOG(0, boost::diagnostic_information(excp));
     530        AllWentWell = false;
     531      }
     532//      LOG(3, "INFO: Orthovector2: " << Orthovector2 << ".")
     533
     534      // create correct coordination for the three atoms
     535      alpha = (TopOrigin->getType()->getHBondAngle(2))/180.*M_PI/2.;  // retrieve triple bond angle from database
     536      l = BondRescale;        // desired bond length
     537      b = 2.*l*sin(alpha);    // base length of isosceles triangle
     538      d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.);   // length for InBondvector
     539      f = b/sqrt(3.);   // length for Orthvector1
     540      g = b/2.;         // length for Orthvector2
     541//      LOG(3, "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", ");
     542//      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));
     543      factors[0] = d;
     544      factors[1] = f;
     545      factors[2] = 0.;
     546      FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
     547      factors[1] = -0.5*f;
     548      factors[2] = g;
     549      SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
     550      factors[2] = -g;
     551      ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
     552
     553      // rescale each to correct BondDistance
     554//      FirstOtherAtom->x.Scale(&BondRescale);
     555//      SecondOtherAtom->x.Scale(&BondRescale);
     556//      ThirdOtherAtom->x.Scale(&BondRescale);
     557
     558      // and relative to *origin atom
     559      *FirstOtherAtom += TopOrigin->getPosition();
     560      *SecondOtherAtom += TopOrigin->getPosition();
     561      *ThirdOtherAtom += TopOrigin->getPosition();
     562
     563      // ... and add to molecule
     564      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
     565      AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
     566      AllWentWell = AllWentWell && AddAtom(ThirdOtherAtom);
     567//      LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
     568//      LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
     569//      LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
     570      Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
     571      Binder->Cyclic = false;
     572      Binder->Type = GraphEdge::TreeEdge;
     573      Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
     574      Binder->Cyclic = false;
     575      Binder->Type = GraphEdge::TreeEdge;
     576      Binder = AddBond(BottomOrigin, ThirdOtherAtom, 1);
     577      Binder->Cyclic = false;
     578      Binder->Type = GraphEdge::TreeEdge;
     579      break;
     580    default:
     581      ELOG(1, "BondDegree does not state single, double or triple bond!");
     582      AllWentWell = false;
     583      break;
     584  }
     585
     586  return AllWentWell;
     587};
    588588
    589589/** Creates a copy of this molecule.

  • src/molecule.hpp

    ree4f2d r2a03b0  
    261261  /// Add/remove atoms to/from molecule.
    262262  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);
    264264  bond::ptr AddBond(atom *first, atom *second, int degree = 1);
    265265  bool hasBondStructure() const;

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

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

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

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

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

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

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

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-anthracene.at

    ree4f2d r2a03b0  
    3535        --select-molecule-by-id 0 \
    3636        --select-molecules-atoms \
    37         --correct-bonddegree \
    3837        --fragment-molecule $FILENAME \
    3938        --order 3 \
    … …  
    4544         --fragment-resultfile ${FILENAME}_results.dat],
    4645         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)/energy) > 1e-5) exit(1)}'], 0)
     46AT_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)
    4847
    4948AT_CLEANUP

  • tests/Calculations/testsuite-calculations-benzene.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-cholesterol.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-coronene.at

    ree4f2d r2a03b0  
    3535        --select-molecule-by-id 0 \
    3636        --select-molecules-atoms \
    37         --correct-bonddegree \
    3837        --fragment-molecule $FILENAME \
    3938        --order 3 \
    … …  
    4544         --fragment-resultfile ${FILENAME}_results.dat],
    4645         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)/energy) > 1e-5) exit(1)}'], 0)
     46AT_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)
    4847
    4948AT_CLEANUP

  • tests/Calculations/testsuite-calculations-cycloheptane.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-dimethyl_bromomalonate.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-glucose.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-heptan.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-isoleucine.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-naphthalene.at

    ree4f2d r2a03b0  
    3535        --select-molecule-by-id 0 \
    3636        --select-molecules-atoms \
    37         --correct-bonddegree \
    3837        --fragment-molecule $FILENAME \
    3938        --order 3 \
    … …  
    4544         --fragment-resultfile ${FILENAME}_results.dat],
    4645         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)/energy) > 1e-5) exit(1)}'], 0)
     46AT_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)
    4847
    4948AT_CLEANUP

  • tests/Calculations/testsuite-calculations-neohexane.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-phenanthrene.at

    ree4f2d r2a03b0  
    3535        --select-molecule-by-id 0 \
    3636        --select-molecules-atoms \
    37         --correct-bonddegree \
    3837        --fragment-molecule $FILENAME \
    3938        --order 3 \
    … …  
    4544         --fragment-resultfile ${FILENAME}_results.dat],
    4645         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)/energy) > 1e-5) exit(1)}'], 0)
     46AT_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)
    4847
    4948AT_CLEANUP

  • tests/Calculations/testsuite-calculations-proline.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-putrescine.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Calculations/testsuite-calculations-tartaric_acid.at

    ree4f2d r2a03b0  
    4545         --fragment-resultfile ${FILENAME}_results.dat],
    4646         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)/energy) > 1e-5) exit(1)}'], 0)
     47AT_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)
    4848
    4949AT_CLEANUP

  • tests/Fragmentations/testsuite.at

    ree4f2d r2a03b0  
    4242
    4343# 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])
    4545
    4646# fragmentation of 1_2-dimethoxyethane

  • tests/Makefile.am

    ree4f2d r2a03b0  
    1313
    1414extracheck:
    15         cd Calculations && $(MAKE) extracheck
     15        cd Calculations && make extracheck
    1616installextracheck:
    17         cd Calculations && $(MAKE) installextracheck
     17        cd Calculations && make installextracheck
Note: See TracChangeset for help on using the changeset viewer.