Changeset 06804b


Ignore:
Timestamp:
Aug 14, 2014, 2:49:51 PM (10 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
3fbdca
Parents:
9fd44f
git-author:
Frederik Heber <heber@…> (08/26/13 19:49:32)
git-committer:
Frederik Heber <heber@…> (08/14/14 14:49:51)
Message:

Removed molecule::AddHydrogenReplacement() and BreadthFirstSearchAdd.

  • this capability has been transfered to SaturatedFragment, and BFSAdd has not been used for a long, long time.
Location:
src
Files:
2 deleted
5 edited

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Exporters/ExportGraph.cpp

    r9fd44f r06804b  
    137137}
    138138
    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 }
     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//}

  • src/Fragmentation/Exporters/ExportGraph.hpp

    r9fd44f r06804b  
    7777  void releaseFragment(SaturatedFragment_ptr &_ptr);
    7878
    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);
     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);
    8383
    8484protected:

  • src/Graph/Makefile.am

    r9fd44f r06804b  
    44GRAPHSOURCE = \
    55        Graph/BondGraph.cpp \
    6         Graph/BreadthFirstSearchAdd.cpp \
    76        Graph/BuildInducedSubgraph.cpp \
    87        Graph/AdjacencyList.cpp \
    … …  
    1312GRAPHHEADER = \
    1413        Graph/BondGraph.hpp \
    15         Graph/BreadthFirstSearchAdd.hpp \
    1614        Graph/BuildInducedSubgraph.hpp \
    1715        Graph/AdjacencyList.hpp \

  • src/molecule.cpp

    r9fd44f r06804b  
    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 };
     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//};
    588588
    589589/** Creates a copy of this molecule.

  • src/molecule.hpp

    r9fd44f r06804b  
    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;
Note: See TracChangeset for help on using the changeset viewer.