Changeset ee4f2d for src/Fragmentation


Ignore:
Timestamp:
Aug 20, 2014, 1:09:40 PM (11 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, Candidate_v1.7.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
fc3aff
Parents:
2a03b0 (diff), 343c5a (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'EnhancingHydrogenSaturation' into stable

Location:
src/Fragmentation
Files:
7 added
12 edited

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Exporters/ExportGraph.cpp

    r2a03b0 ree4f2d  
    5858  const Graph &_graph,
    5959  const enum HydrogenTreatment _treatment,
    60   const enum HydrogenSaturation _saturation) :
     60  const enum HydrogenSaturation _saturation,
     61  const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
    6162  TotalGraph(_graph),
    6263  BondFragments(World::getPointer()),
    6364  treatment(_treatment),
    6465  saturation(_saturation),
     66  globalsaturationpositions(_globalsaturationpositions),
    6567  CurrentKeySet(TotalGraph.begin())
    6668{
    … …  
    115117          hydrogens,
    116118          treatment,
    117           saturation)
     119          saturation,
     120          globalsaturationpositions)
    118121  );
    119122  // and return
    … …  
    137140}
    138141
    139 /** Internal helper to create from each keyset a molecule
    140  *
    141  */
    142 void 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  */
    163 molecule * 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  */
    184 int 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  */
    204 void 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 }
     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//}

  • src/Fragmentation/Exporters/ExportGraph.hpp

    r2a03b0 ree4f2d  
    4141      const Graph &_graph,
    4242      const enum HydrogenTreatment _treatment,
    43       const enum HydrogenSaturation _saturation);
     43      const enum HydrogenSaturation _saturation,
     44      const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
    4445  virtual ~ExportGraph();
    4546
    … …  
    7778  void releaseFragment(SaturatedFragment_ptr &_ptr);
    7879
    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);
     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);
    8384
    8485protected:
    … …  
    102103  const enum HydrogenSaturation saturation;
    103104
     105  //!> Global information over all atoms with saturation positions to be used per fragment
     106  const SaturatedFragment::GlobalSaturationPositions_t &globalsaturationpositions;
     107
    104108private:
    105109  //!> iterator pointing at the CurrentKeySet to be exported

  • src/Fragmentation/Exporters/ExportGraph_ToFiles.cpp

    r2a03b0 ree4f2d  
    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
    5961 */
    6062ExportGraph_ToFiles::ExportGraph_ToFiles(
    6163    const Graph &_graph,
    6264    const enum HydrogenTreatment _treatment,
    63     const enum HydrogenSaturation _saturation) :
    64                 ExportGraph(_graph, _treatment, _saturation)
     65    const enum HydrogenSaturation _saturation,
     66    const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
     67                ExportGraph(_graph, _treatment, _saturation, _globalsaturationpositions)
    6568{}
    6669

  • src/Fragmentation/Exporters/ExportGraph_ToFiles.hpp

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

  • src/Fragmentation/Exporters/ExportGraph_ToJobs.cpp

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

  • src/Fragmentation/Exporters/ExportGraph_ToJobs.hpp

    r2a03b0 ree4f2d  
    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
    3739   */
    3840  ExportGraph_ToJobs(
    3941      const Graph &_graph,
    4042      const enum HydrogenTreatment _treatment,
    41       const enum HydrogenSaturation _saturation);
     43      const enum HydrogenSaturation _saturation,
     44      const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions);
    4245  virtual ~ExportGraph_ToJobs();
    4346

  • src/Fragmentation/Exporters/SaturatedFragment.cpp

    r2a03b0 ree4f2d  
    6363    HydrogenPool &_hydrogens,
    6464    const enum HydrogenTreatment _treatment,
    65     const enum HydrogenSaturation _saturation) :
     65    const enum HydrogenSaturation _saturation,
     66    const GlobalSaturationPositions_t &_globalsaturationpositions) :
    6667  container(_container),
    6768  set(_set),
    … …  
    7778  container.insert(set);
    7879
    79   // prepare saturation hydrogens
    80   saturate();
     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);
    8186}
    8287
    … …  
    99104}
    100105
    101 void 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();
     106typedef std::vector<atom *> atoms_t;
     107
     108atoms_t gatherAllAtoms(const KeySet &_FullMolecule)
     109{
     110  atoms_t atoms;
     111  for (KeySet::const_iterator iter = _FullMolecule.begin();
     112      iter != _FullMolecule.end();
    107113      ++iter) {
    108114    atom * const Walker = World::getInstance().getAtom(AtomById(*iter));
    109115    ASSERT( Walker != NULL,
    110         "SaturatedFragment::OutputConfig() - id "
     116        "gatherAllAtoms() - id "
    111117        +toString(*iter)+" is unknown to World.");
    112118    atoms.push_back(Walker);
    113119  }
    114120
    115 //  bool LonelyFlag = false;
    116   for (std::vector<atom *>::const_iterator iter = atoms.begin();
    117       iter != atoms.end();
     121  return atoms;
     122}
     123
     124typedef std::map<atom *, BondList > CutBonds_t;
     125
     126CutBonds_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();
    118135      ++iter) {
    119136    atom * const Walker = *iter;
    … …  
    125142        ++BondRunner) {
    126143      atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
    127       // if in set
    128       if (set.find(OtherWalker->getId()) != set.end()) {
     144      // if other atom is in key set or is a specially treated hydrogen
     145      if (_set.find(OtherWalker->getId()) != _set.end()) {
    129146        LOG(4, "DEBUG: Walker " << *Walker << " is bound to " << *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;
     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 << ".");
    140151      } else {
    141152        LOG(4, "DEBUG: Walker " << *Walker << " is bound to "
    142153            << *OtherWalker << ", who is not in this fragment molecule.");
    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);
     154          if (CutBonds.count(Walker) == 0)
     155            CutBonds.insert( std::make_pair(Walker, BondList() ));
     156          CutBonds[Walker].push_back(*BondRunner);
     157      }
     158    }
     159  }
     160
     161  return CutBonds;
     162}
     163
     164typedef std::vector<atomId_t> atomids_t;
     165
     166atomids_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
     200void 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
     224void 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());
    151295        }
    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();
    158296      }
    159297    }
    160   }
     298  } else
     299    LOG(3, "INFO: We are not saturating cut bonds.");
     300}
     301
     302const 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
     319bool 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;
    161335}
    162336
    … …  
    270444  if (BondRescale == -1) {
    271445    ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->getDegree() << "!");
    272     return false;
    273     BondRescale = bondlength;
     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    }
    274452  } else {
    275453    if (!IsAngstroem)

  • src/Fragmentation/Exporters/SaturatedFragment.hpp

    r2a03b0 ree4f2d  
    1818#include <string>
    1919
     20#include "Atom/atom_bondedparticleinfo.hpp"
    2021#include "Bond/bond.hpp"
    2122#include "Fragmentation/KeySet.hpp"
    2223#include "Fragmentation/HydrogenSaturation_enum.hpp"
    2324#include "Parser/FormatParserStorage.hpp"
     25
     26#include "LinearAlgebra/Vector.hpp"
    2427
    2528class atom;
    … …  
    4245  typedef std::set<KeySet> KeySetsInUse_t;
    4346
     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
    4454  /** Constructor of SaturatedFragment requires \a set which we are tightly
    4555   * associated.
    … …  
    4858   * \param _container container to add KeySet as in-use
    4959   * \param _hydrogens pool with hydrogens for saturation
     60   * \param _globalsaturationpositions saturation positions to be used
    5061   */
    5162  SaturatedFragment(
    … …  
    5465      HydrogenPool &_hydrogens,
    5566      const enum HydrogenTreatment _treatment,
    56       const enum HydrogenSaturation saturation);
     67      const enum HydrogenSaturation saturation,
     68      const GlobalSaturationPositions_t &_globalsaturationpositions);
    5769
    5870  /** Destructor of class SaturatedFragment.
    … …  
    100112  /** Helper function to lease and bring in place saturation hydrogens.
    101113   *
     114   * Here, we use local information to calculate saturation positions.
     115   *
    102116   */
    103117  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);
    104135
    105136  /** Helper function to get a hydrogen replacement for a given \a replacement.
    … …  
    109140   */
    110141  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);
    111155
    112156  /** Leases and adds a Hydrogen atom in replacement for the given atom \a *partner in bond with a *origin.

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

    r2a03b0 ree4f2d  
    44FRAGMENTATIONEXPORTERSSOURCES = \
    55        ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.cpp \
    6         ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp
     6        ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.cpp \
     7        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.cpp
    78
    89FRAGMENTATIONEXPORTERSTESTSHEADERS = \
    910        ../Fragmentation/Exporters/unittests/HydrogenPoolUnitTest.hpp \
    10         ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.hpp
     11        ../Fragmentation/Exporters/unittests/SaturatedFragmentUnitTest.hpp \
     12        ../Fragmentation/Exporters/unittests/SaturationDistanceMaximizerUnitTest.hpp
    1113
    1214FRAGMENTATIONEXPORTERSTESTS = \
    1315        HydrogenPoolUnitTest \
    14         SaturatedFragmentUnitTest
     16        SaturatedFragmentUnitTest \
     17        SaturationDistanceMaximizerUnitTest
    1518
    1619TESTS += $(FRAGMENTATIONEXPORTERSTESTS)
    … …  
    3740        $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
    3841        ${FRAGMENTATIONEXPORTERSLIBS}
     42
     43SaturationDistanceMaximizerUnitTest_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
     47SaturationDistanceMaximizerUnitTest_LDADD = \
     48        ../libMolecuilderUI.la  \
     49        $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
     50        ${FRAGMENTATIONEXPORTERSLIBS}
    3951       
    4052

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

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

  • src/Fragmentation/Makefile.am

    r2a03b0 ree4f2d  
    77        Fragmentation/Exporters/ExportGraph.cpp \
    88        Fragmentation/Exporters/HydrogenPool.cpp \
     9        Fragmentation/Exporters/SaturatedBond.cpp \
    910        Fragmentation/Exporters/SaturatedFragment.cpp \
     11        Fragmentation/Exporters/SaturationDistanceMaximizer.cpp \
    1012        Fragmentation/Homology/FragmentEdge.cpp \
    1113        Fragmentation/Homology/FragmentNode.cpp \
    … …  
    3335        Fragmentation/Exporters/ExportGraph.hpp \
    3436        Fragmentation/Exporters/HydrogenPool.hpp \
     37        Fragmentation/Exporters/SaturatedBond.hpp \
    3538        Fragmentation/Exporters/SaturatedFragment.hpp \
     39        Fragmentation/Exporters/SaturationDistanceMaximizer.hpp \
    3640        Fragmentation/Homology/FragmentEdge.hpp \
    3741        Fragmentation/Homology/FragmentNode.hpp \

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

    r2a03b0 ree4f2d  
    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    }
    115123    IndexedVectors::indices_t indices(arrayiter->begin(), arrayiter->end());
    116124    boost::fusion::at_key<MPQCDataFused::forces>(instance) =
Note: See TracChangeset for help on using the changeset viewer.