Changeset 266237

Timestamp:

Oct 18, 2009, 2:51:38 PM (15 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:

174e0e

Parents:

b8b75d

git-author:

Frederik Heber <heber@…> (10/18/09 14:15:37)

git-committer:

Frederik Heber <heber@…> (10/18/09 14:51:38)

Message:

Huge refactoring: molecule::ListOfBondsPerAtom and molecule::NumberOfBondsPerAtom removed, atom::ListOfBonds introduced. Unit Test for ListOfBonds manipulation introduced.

• changes to class atom:
  • ~atom(): removes all bonds it is part of and unlink()s
  • new functions: CountBonds() - counts bonds times BondDegree, OutputGraphInfo() - graph data is output, OutputComponentNumber() - output of component numbers, OutputAdjacency() - indices of this and its neighbours, CorrectBondDegree() - corrects bond degree incrementally (by at most one) -- all since ListOfBonds is now local info
  • OutputBondOfAtom() - does not need ListOfBondsPerAtom, ... as ListOfBonds is now local info
  • new functions for ListOfBonds access: RegisterBond(), UnregisterBond(), UnregisterAllBonds() (the latter only remove from the list, do not delete, as molecule still has a global list of all bonds and does deletion.)

• changes to class bond:
  • ~bond() uses unlink(), uses UnregisterBond() on leftatom and rightatom
  • removed unused bond::GetFirstBond() and bond::GetLastBond()
  • GetOtherAtom() has now const parameter

• changes to builder.cpp: removed CreateListOfBondsPerAtom() calls, as the creation of the global arrays is not necessary anymore

• changes to LinkedCell: LinkedCell::CheckBounds(int[NDIM]) does not admonish out of bonds as this is not desired for the local offset which may become out of bounds.

• changes to lists.hpp templates: BUGFIX: unlink() now sets ->next and ->previous to NULL, cleanup() uses removedwithoutcheck()

• new templates for molecule.hpp: SumPerAtom() allows for summation of the return value of atom:...() member fiunctions. This is needed e.g. for atom::CorrectBondDegree()

• changes to MoleculeListClass:
  • changes due to ListOfBondsPerAtom -> atom::ListOfBonds: AddHydrogenCorrection(), FillBondStructureFromReference()

• changes to class molecule:
  • changes due to ListOfBondsPerAtom -> atom::ListOfBonds: GuesstimateFragmentCount(), FragmentMolecule(), FillSPListandLabelVertices(), ScanForPeriodicCorrection(), DeterminePeriodicCenter(), CyclicStructureAnalysis(), InitComponentNumbers(), FindNextUnused(), CheckAdjacencyFileAgainstMolecule(), PickLocalBackEdges(), BreadthFirstSearchAdd(), BuildInducedSubgraph(), CheckForConnectedSubgraph(), molecule(), ~molecule(), AddHydrogenReplacementAtom(), RemoveBonds(), OutputListOfBonds(), removed CreateListOfBondsPerAtom()
  • due to removal of CreateListOfBonds(): FragmentMolecule(), CreateAdjacencyListFromDbondFile(), CreateAdjacencyList(), SetNextComponentNumber()
  • ActOnAllAtoms() used: CreateAdjacencyList() - uses atom::OutputBondOfAtom(), StoreAdjacencyToFile() - uses atom::OutputAdjacency(),
  • SumPerAtom() used: CorrectBondDegree() - uses atom::CorrectBondDegree(),
  • function removed: molecule::CountAtomsBonds() - shifted to atom::CountBonds(), OutputComponentNumber() - shifted to atom::OutputComponentNumber()
  • CountCyclicBonds() - manipulates molecule::NoCyclicBonds directly.
  • rewritten DepthFirstSearchAnalysis() - refactored to CyclicBondAnalysis(), OutputGraphInfoPerAtom(), OutputGraphInfoPerBond()
  • rewritten RemoveAtom() - uses RemoveBonds()
  • FIX: rewritten CleanupMolecule() - cleanup was wrong order, now is first bonds, then atoms
  • rewritten BreadthFirstSearchAdd() - uses molecule::CopyBond() (equal code was present multiple(!) times before)
  • new functions CyclicBondAnalysis(), OutputGraphInfoPerAtom(), OutputGraphInfoPerBond(), CopyBond()

• New UnitTest for ListOfBonds:
  • tests AddBond(), RemoveAtom(), RemoveBond(), RemoveBonds() and simple deletion of atoms or bonds, where afterwards chain lists of bonds should remain intact.

Signed-off-by: Frederik Heber <heber@…>

Location:

src

Files:

• atom.cpp (modified) (7 diffs)
• atom.hpp (modified) (5 diffs)
• bond.cpp (modified) (4 diffs)
• bond.hpp (modified) (2 diffs)
• builder.cpp (modified) (2 diffs)
• linkedcell.cpp (modified) (2 diffs)
• lists.hpp (modified) (2 diffs)
• molecule.cpp (modified) (12 diffs)
• molecule.hpp (modified) (6 diffs)
• molecule_fragmentation.cpp (modified) (14 diffs)
• molecule_geometry.cpp (modified) (3 diffs)
• molecule_graph.cpp (modified) (30 diffs)
• molecule_template.hpp (modified) (1 diff)
• moleculelist.cpp (modified) (5 diffs)
• unittests/Makefile.am (modified) (2 diffs)
• unittests/listofbondsunittest.cpp (added)
• unittests/listofbondsunittest.hpp (added)

src/atom.cpp

@@ -9 +9 @@
 #include "config.hpp"
 #include "element.hpp"
+#include "lists.hpp"
 #include "memoryallocator.hpp"
 #include "parser.hpp"
@@ -67 +68 @@
 atom::~atom()
 {
+  BondList::const_iterator Runner;
+  while (!ListOfBonds.empty()) {
+    Runner = ListOfBonds.begin();
+    removewithoutcheck(*Runner);
+  }
+  unlink(this);
   Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
   Free<char>(&Name, "atom::~atom: *Name");
@@ -118 +125 @@
 {
   return (node->IsInParallelepiped(offset, parallelepiped));
+};
+
+/** Counts the number of bonds weighted by bond::BondDegree.
+ * \param bonds times bond::BondDegree
+ */
+int atom::CountBonds() const
+{
+  int NoBonds = 0;
+  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
+    NoBonds += (*Runner)->BondDegree;
+  return NoBonds;
 };
 
@@ -215 +233 @@
 };
 
-/** Prints all bonds of this atom from given global lists.
+/** Output graph info of this atom.
+ * \param *out output stream
+ */
+void atom::OutputGraphInfo(ofstream *out) const
+{
+  *out << Verbose(2) << "Atom " << Name << " is " << ((SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
+  OutputComponentNumber(out);
+  *out << " with Lowpoint " << LowpointNr << " and Graph Nr. " << GraphNr << "." << endl;
+};
+
+/** Output a list of flags, stating whether the bond was visited or not.
+ * Note, we make use of the last entry of the ComponentNr always being -1 if allocated.
+ * \param *out output stream for debugging
+ */
+void atom::OutputComponentNumber(ofstream *out) const
+{
+  if (ComponentNr != NULL) {
+    for (int i=0; ComponentNr[i] != -1; i++)
+      *out << ComponentNr[i] << " ";
+  }
+};
+
+
+/** Prints all bonds of this atom with total degree.
  * \param *out stream to output to
- * \param *NumberOfBondsPerAtom array with number of bonds per atomic index
- * \param ***ListOfBondsPerAtom array per atomic index of array with pointer to bond
  * \return true - \a *out present, false - \a *out is NULL
  */
-bool atom::OutputBondOfAtom(ofstream *out, int *NumberOfBondsPerAtom, bond ***ListOfBondsPerAtom) const
+bool atom::OutputBondOfAtom(ofstream *out) const
 {
   if (out != NULL) {
@@ -227 +266 @@
     if (type->Z != 1) {   // regard only non-hydrogen
 #endif
-      *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << NumberOfBondsPerAtom[nr] << " bonds: ";
+      *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << ListOfBonds.size() << " bonds: ";
       int TotalDegree = 0;
-      for (int j=0;j<NumberOfBondsPerAtom[nr];j++) {
-        *out << *ListOfBondsPerAtom[nr][j] << "\t";
-        TotalDegree += ListOfBondsPerAtom[nr][j]->BondDegree;
+      for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) {
+        *out << **Runner << "\t";
+        TotalDegree += (*Runner)->BondDegree;
       }
       *out << " -- TotalDegree: " << TotalDegree << endl;
@@ -240 +279 @@
   } else
     return false;
+};
+
+/** Output of atom::nr along with all bond partners.
+ * \param *AdjacencyFile output stream
+ */
+void atom::OutputAdjacency(ofstream *AdjacencyFile) const
+{
+  *AdjacencyFile << nr << "\t";
+  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
+    *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
+  *AdjacencyFile << endl;
+};
+
+/** Puts a given bond into atom::ListOfBonds.
+ * \param *Binder bond to insert
+ */
+bool atom::RegisterBond(bond *Binder)
+{
+  bool status = false;
+  if (Binder != NULL) {
+    if (Binder->Contains(this)) {
+      ListOfBonds.push_back(Binder);
+      status = true;
+    } else {
+      cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
+    }
+  } else {
+    cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
+  }
+  return status;
+};
+
+/** Removes a given bond from atom::ListOfBonds.
+ * \param *Binder bond to remove
+ */
+bool atom::UnregisterBond(bond *Binder)
+{
+  bool status = false;
+  if (Binder != NULL) {
+    if (Binder->Contains(this)) {
+      ListOfBonds.remove(Binder);
+      status = true;
+    } else {
+      cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
+    }
+  } else {
+    cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
+  }
+  return status;
+};
+
+/** Removes all bonds from atom::ListOfBonds.
+ * \note Does not do any memory de-allocation.
+ */
+void atom::UnregisterAllBond()
+{
+  ListOfBonds.clear();
+};
+
+/** Corrects the bond degree by one at most if necessary.
+ * \param *out output stream for debugging
+ */
+int atom::CorrectBondDegree(ofstream *out)
+{
+  int NoBonds = 0;
+  int OtherNoBonds = 0;
+  int FalseBondDegree = 0;
+  atom *OtherWalker = NULL;
+  bond *CandidateBond = NULL;
+
+  *out << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+  NoBonds = CountBonds();
+  if ((int)(type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
+    for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
+      OtherWalker = (*Runner)->GetOtherAtom(this);
+      OtherNoBonds = OtherWalker->CountBonds();
+      *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+      if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
+        if ((CandidateBond == NULL) || (ListOfBonds.size() > OtherWalker->ListOfBonds.size())) { // pick the one with fewer number of bonds first
+          CandidateBond = (*Runner);
+          *out << Verbose(3) << "New candidate is " << *CandidateBond << "." << endl;
+        }
+      }
+    }
+    if ((CandidateBond != NULL)) {
+      CandidateBond->BondDegree++;
+      *out << Verbose(2) << "Increased bond degree for bond " << *CandidateBond << "." << endl;
+    } else {
+      *out << Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl;
+      FalseBondDegree++;
+    }
+  }
+  return FalseBondDegree;
 };
 
@@ -345 +477 @@
 {
   *file << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << endl;
-  //cout << Verbose(2) << "Storing: " << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "\t" << (int)Walker->MaxOrder << "." << endl;
+  //cout << Verbose(2) << "Storing: " << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << "." << endl;
 }
 

src/atom.hpp

@@ -19 +19 @@
 
 #include <iostream>
+#include <list>
 #include <vector>
 
@@ -32 +33 @@
 class ForceMatrix;
 class Vector;
+
+#define BondList list<bond *>
 
 /********************************************** declarations *******************************/
@@ -50 +53 @@
     Vector v;       //!< velocity vector of atom, giving last velocity within cell
     Vector F;       //!< Force vector of atom, giving last force within cell
+    BondList ListOfBonds; //!< list of all bonds
     element *type;  //!< pointing to element
     atom *previous; //!< previous atom in molecule list
@@ -76 +80 @@
   bool OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const;
   bool OutputTrajectoryXYZ(ofstream *out, int step) const;
-  bool OutputBondOfAtom(ofstream *out, int *NumberOfBondsPerAtom, bond ***ListOfBondsPerAtom) const;
+  bool OutputBondOfAtom(ofstream *out) const;
+  void OutputAdjacency(ofstream *AdjacencyFile) const;
+
 
   void EqualsFather ( atom *ptr, atom **res );
@@ -96 +102 @@
   // bond order stuff
   void OutputOrder(ofstream *file);
+  void OutputGraphInfo(ofstream *out) const;
+  void OutputComponentNumber(ofstream *out) const;
+  int CountBonds() const;
+  int CorrectBondDegree(ofstream *out);
+
+  bool RegisterBond(bond *Binder);
+  bool UnregisterBond(bond *Binder);
+  void UnregisterAllBond();
 
   // constraint potential and dynamics stuff

src/bond.cpp

@@ -53 +53 @@
   Cyclic = false;
 };
+
 bond::bond(atom *left, atom *right)
 {
@@ -77 +78 @@
 {
   // remove this node from the list structure
-  if (previous != NULL) {
-    previous->next = next;
-  }
-  if (next != NULL) {
-    next->previous = previous;
-  }
+  if (leftatom != NULL)
+    leftatom->UnregisterBond(this);
+  if (rightatom != NULL)
+  rightatom->UnregisterBond(this);
+  unlink(this);
 };
 
@@ -95 +95 @@
  * \return pointer to the other atom in the bond, NULL if no match (indicates something's wrong with the bond) 
  */
-atom * bond::GetOtherAtom(atom *Atom) const
+atom * bond::GetOtherAtom(const atom *Atom) const
 {
   if(leftatom == Atom) 
@@ -105 +105 @@
 };
 
-/** Get the other atom in a bond if one is specified.
- * \param *Atom the pointer to the one atom
- * \return pointer to the other atom in the bond, NULL if no match (indicates something's wrong with the bond) 
- */
-bond * bond::GetFirstBond()
-{
-  return GetFirst(this);
-};
-
-/** Get the other atom in a bond if one is specified.
- * \param *Atom the pointer to the one atom
- * \return pointer to the other atom in the bond, NULL if no match (indicates something's wrong with the bond) 
- */
-bond * bond::GetLastBond()
-{
-  return GetLast(this);
-};
 
 /** Returns whether vertex was used in DFS.

src/bond.hpp

@@ -25 +25 @@
 
 /** Bonds between atoms.
- * Class incorporates bonds between atoms in a molecule,
- * used to derive tge fragments in many-body bond order
- * calculations.
+ * Class incorporates bonds between atoms in a molecule.
+ * Note that we regard bond always as something in a molecule,
+ * as it is the glue making up the connected subgrapgh and
+ * hence the molecule. Thus, bonds belong globally to the molecule
+ * (and are free'd there) and only locally to the atom classs.
  */
 class bond {
@@ -41 +43 @@
     enum EdgeType Type;//!< whether this is a tree or back edge
 
-  atom * GetOtherAtom(atom *Atom) const;
-  bond * GetFirstBond();
-  bond * GetLastBond();
+  atom * GetOtherAtom(const atom *Atom) const;
 
   bool MarkUsed(enum Shading color);

src/builder.cpp

@@ -895 +895 @@
           start = clock();
           mol->CreateAdjacencyList((ofstream *)&cout, bonddistance, configuration->GetIsAngstroem());
-          mol->CreateListOfBondsPerAtom((ofstream *)&cout);
           end = clock();
           cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
@@ -1555 +1554 @@
                 class StackClass<bond *> *LocalBackEdgeStack = NULL;
                 mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem());
-                mol->CreateListOfBondsPerAtom((ofstream *)&cout);
                 Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
                 if (Subgraphs != NULL) {

src/linkedcell.cpp

@@ -203 +203 @@
 
 /** Checks whether LinkedCell::n[] plus relative offset is each within [0,N[]].
+ * Note that for this check we don't admonish if out of bounds.
  * \param relative[NDIM] relative offset to current cell
  * \return if all in intervals - true, else -false
@@ -211 +212 @@
   for(int i=0;i<NDIM;i++)
     status = status && ((n[i]+relative[i] >=0) && (n[i]+relative[i] < N[i]));
-  if (!status)
-  cerr << "ERROR: indices are out of bounds!" << endl;
   return status;
 };

src/lists.hpp

@@ -35 +35 @@
   if (walker->previous != NULL)
     walker->previous->next = walker->next;
+  walker->next = NULL;
+  walker->previous= NULL;
 };
 
@@ -112 +114 @@
 {
   X *pointer = start->next;
-  X *walker;
+  X *walker = NULL;
   while (pointer != end) { // go through list
     walker = pointer; // mark current
     pointer = pointer->next; // step onward beforehand
     // remove walker
-    unlink(walker);
-    delete(walker);
-    walker = NULL;
+    removewithoutcheck(walker);
   }
   return true;

src/molecule.cpp

@@ -48 +48 @@
   NoNonHydrogen = 0;
   NoCyclicBonds = 0;
-  ListOfBondsPerAtom = NULL;
-  NumberOfBondsPerAtom = NULL;
   ElementCount = 0;
   for(int i=MAX_ELEMENTS;i--;)
@@ -66 +64 @@
 molecule::~molecule()
 {
-  if (ListOfBondsPerAtom != NULL)
-    for(int i=AtomCount;i--;)
-      Free(&ListOfBondsPerAtom[i]);
-  Free(&ListOfBondsPerAtom);
-  Free(&NumberOfBondsPerAtom);
   if (TesselStruct != NULL)
     delete(TesselStruct);
@@ -160 +153 @@
  * \param *origin pointer to atom which acts as the origin for scaling the added hydrogen to correct bond length
  * \param *replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
- * \param **BondList list of bonds \a *replacement has (necessary to determine plane for double and triple bonds)
- * \param NumBond  number of bonds in \a **BondList
  * \param isAngstroem whether the coordination of the given atoms is in AtomicLength (false) or Angstrom(true)
  * \return number of atoms added, if < bond::BondDegree then something went wrong
  * \todo double and triple bonds splitting (always use the tetraeder angle!)
  */
-bool molecule::AddHydrogenReplacementAtom(ofstream *out, bond *TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem)
+bool molecule::AddHydrogenReplacementAtom(ofstream *out, bond *TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
 {
   double bondlength;  // bond length of the bond to be replaced/cut
@@ -177 +168 @@
   Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
   Vector InBondvector;    // vector in direction of *Bond
+  double *matrix;
   bond *Binder = NULL;
-  double *matrix;
 
 //  *out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
@@ -248 +239 @@
     case 2:
       // determine two other bonds (warning if there are more than two other) plus valence sanity check
-      for (int i=0;i<NumBond;i++) {
-        if (BondList[i] != TopBond) {
+      for (BondList::const_iterator Runner = TopOrigin->ListOfBonds.begin(); Runner != TopOrigin->ListOfBonds.end(); (++Runner)) {
+        if ((*Runner) != TopBond) {
           if (FirstBond == NULL) {
-            FirstBond = BondList[i];
-            FirstOtherAtom = BondList[i]->GetOtherAtom(TopOrigin);
+            FirstBond = (*Runner);
+            FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
           } else if (SecondBond == NULL) {
-            SecondBond = BondList[i];
-            SecondOtherAtom = BondList[i]->GetOtherAtom(TopOrigin);
+            SecondBond = (*Runner);
+            SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
           } else {
             *out << Verbose(3) << "WARNING: Detected more than four bonds for atom " << TopOrigin->Name;
@@ -523 +514 @@
   if (first->next != last) {  // if adjaceny list is present
     copy->BondDistance = BondDistance;
-    copy->CreateListOfBondsPerAtom((ofstream *)&cout);
   }
 
@@ -557 +547 @@
   if ((atom1 != NULL) && (FindAtom(atom1->nr) != NULL) && (atom2 != NULL) && (FindAtom(atom2->nr) != NULL)) {
     Binder = new bond(atom1, atom2, degree, BondCount++);
+    atom1->RegisterBond(Binder);
+    atom2->RegisterBond(Binder);
     if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
       NoNonBonds++;
@@ -585 +577 @@
 bool molecule::RemoveBonds(atom *BondPartner)
 {
-  cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
+  //cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
+  BondList::const_iterator ForeRunner;
+  while (!BondPartner->ListOfBonds.empty()) {
+    ForeRunner = BondPartner->ListOfBonds.begin();
+    RemoveBond(*ForeRunner);
+  }
   return false;
 };
@@ -635 +632 @@
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
+  RemoveBonds(pointer);
   return remove(pointer, start, end);
 };
@@ -661 +659 @@
 bool molecule::CleanupMolecule()
 {
-  return (cleanup(start,end) && cleanup(first,last));
+  return (cleanup(first,last) && cleanup(start,end));
 };
 
@@ -773 +771 @@
 };
 
-/** Outputs contents of molecule::ListOfBondsPerAtom.
+/** Outputs contents of each atom::ListOfBonds.
  * \param *out output stream
  */
 void molecule::OutputListOfBonds(ofstream *out) const
 {
-  *out << Verbose(2) << endl << "From Contents of ListOfBondsPerAtom, all non-hydrogen atoms:" << endl;
-  ActOnAllAtoms (&atom::OutputBondOfAtom, out, NumberOfBondsPerAtom, ListOfBondsPerAtom);
+  *out << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl;
+  ActOnAllAtoms (&atom::OutputBondOfAtom, out);
   *out << endl;
 };
@@ -900 +898 @@
   }
   configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.MaxPsiDouble;
-};
-
-
-/** Creates an 2d array of pointer with an entry for each atom and each bond it has.
- * Updates molecule::ListOfBondsPerAtom, molecule::NumberOfBondsPerAtom by parsing through
- * bond chain list, using molecule::AtomCount and molecule::BondCount.
- * Allocates memory, fills the array and exits
- * \param *out output stream for debugging
- */
-void molecule::CreateListOfBondsPerAtom(ofstream *out)
-{
-  bond *Binder = NULL;
-  *out << Verbose(1) << "Begin of Creating ListOfBondsPerAtom: AtomCount = " << AtomCount << "\tBondCount = " << BondCount << "\tNoNonBonds = " << NoNonBonds << "." << endl;
-
-  // re-allocate memory
-  *out << Verbose(2) << "(Re-)Allocating memory." << endl;
-  if (ListOfBondsPerAtom != NULL) {
-    for(int i=AtomCount;i--;)
-      Free(&ListOfBondsPerAtom[i]);
-    Free(&ListOfBondsPerAtom);
-  }
-  if (NumberOfBondsPerAtom != NULL)
-    Free(&NumberOfBondsPerAtom);
-  ListOfBondsPerAtom = Malloc<bond**>(AtomCount, "molecule::CreateListOfBondsPerAtom: ***ListOfBondsPerAtom");
-  NumberOfBondsPerAtom = Malloc<int>(AtomCount, "molecule::CreateListOfBondsPerAtom: *NumberOfBondsPerAtom");
-
-  // reset bond counts per atom
-  for(int i=AtomCount;i--;)
-    NumberOfBondsPerAtom[i] = 0;
-  // count bonds per atom
-  Binder = first;
-  while (Binder->next != last) {
-    Binder = Binder->next;
-    NumberOfBondsPerAtom[Binder->leftatom->nr]++;
-    NumberOfBondsPerAtom[Binder->rightatom->nr]++;
-  }
-  for(int i=AtomCount;i--;) {
-    // allocate list of bonds per atom
-    ListOfBondsPerAtom[i] = Malloc<bond*>(NumberOfBondsPerAtom[i], "molecule::CreateListOfBondsPerAtom: **ListOfBondsPerAtom[]");
-    // clear the list again, now each NumberOfBondsPerAtom marks current free field
-    NumberOfBondsPerAtom[i] = 0;
-  }
-  // fill the list
-  Binder = first;
-  while (Binder->next != last) {
-    Binder = Binder->next;
-    ListOfBondsPerAtom[Binder->leftatom->nr][NumberOfBondsPerAtom[Binder->leftatom->nr]++] = Binder;
-    ListOfBondsPerAtom[Binder->rightatom->nr][NumberOfBondsPerAtom[Binder->rightatom->nr]++] = Binder;
-  }
-
-  // output list for debugging
-  *out << Verbose(3) << "ListOfBondsPerAtom for each atom:" << endl;
-  ActOnAllAtoms( &atom::OutputBondOfAtom, out, NumberOfBondsPerAtom, ListOfBondsPerAtom );
-
-  *out << Verbose(1) << "End of Creating ListOfBondsPerAtom." << endl << endl;
 };
 

src/molecule.hpp

@@ -86 +86 @@
     bond *first;        //!< start of bond list
     bond *last;         //!< end of bond list
-    bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
     int MDSteps;        //!< The number of MD steps in Trajectories
-    int *NumberOfBondsPerAtom;  //!< Number of Bonds each atom has
     int AtomCount;          //!< number of atoms, brought up-to-date by CountAtoms()
     int BondCount;          //!< number of atoms, brought up-to-date by CountBonds()
@@ -190 +188 @@
   template <typename T> void SetAtomValueToIndexedArray ( T *array, int element::*index, T atom::*value );
 
+  template <typename res> res SumPerAtom(res (atom::*f)() );
+  template <typename res> res SumPerAtom(res (atom::*f)() const );
+  template <typename res> res SumPerAtom(res (atom::*f)() ) const;
+  template <typename res> res SumPerAtom(res (atom::*f)() const ) const;
+  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) , T t );
+  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) const, T t );
+  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) , T t ) const;
+  template <typename T, typename res> res SumPerAtom(res (atom::*f)(T) const, T t ) const;
+
   /// remove atoms from molecule.
   bool AddAtom(atom *pointer);
@@ -199 +206 @@
   atom * AddCopyAtom(atom *pointer);
   bool AddXYZFile(string filename);
-  bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);
+  bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem);
   bond * AddBond(atom *first, atom *second, int degree = 1);
   bool RemoveBond(bond *pointer);
@@ -247 +254 @@
   void CreateAdjacencyListFromDbondFile(ofstream *out, ifstream *output);
   void CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem);
-  void CreateListOfBondsPerAtom(ofstream *out);
   int CorrectBondDegree(ofstream *out);
   void OutputBondsList(ofstream *out);
-  int CountAtomsBonds(int nr);
+  void CyclicBondAnalysis();
+  void OutputGraphInfoPerAtom(ofstream *out);
+  void OutputGraphInfoPerBond(ofstream *out);
 
 
@@ -260 +268 @@
   void SetNextComponentNumber(atom *vertex, int nr);
   void InitComponentNumbers();
-  void OutputComponentNumber(ofstream *out, atom *vertex);
   void ResetAllBondsToUnused();
   void ResetAllAtomNumbers();
@@ -266 +273 @@
   bool CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment);
   string GetColor(enum Shading color);
+  bond * CopyBond(atom *left, atom *right, bond *CopyBond);
+
 
   molecule *CopyMolecule();

src/molecule_fragmentation.cpp

@@ -31 +31 @@
 int molecule::GuesstimateFragmentCount(ofstream *out, int order)
 {
-  int c = 0;
+  size_t c = 0;
   int FragmentCount;
   // get maximum bond degree
@@ -37 +37 @@
   while (Walker->next != end) {
     Walker = Walker->next;
-    c = (NumberOfBondsPerAtom[Walker->nr] > c) ? NumberOfBondsPerAtom[Walker->nr] : c;
+    c = (Walker->ListOfBonds.size() > c) ? Walker->ListOfBonds.size() : c;
   }
   FragmentCount = NoNonHydrogen*(1 << (c*order));
@@ -564 +564 @@
 
   // ===== 1. Check whether bond structure is same as stored in files ====
-
-  // fill the adjacency list
-  CreateListOfBondsPerAtom(out);
 
   // create lookup table for Atom::nr
@@ -630 +627 @@
       MolecularWalker = MolecularWalker->next;
       *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
-      //MolecularWalker->Leaf->OutputListOfBonds(out);  // output ListOfBondsPerAtom for debugging
+      //MolecularWalker->Leaf->OutputListOfBonds(out);  // output atom::ListOfBonds for debugging
       if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
         // call BOSSANOVA method
@@ -879 +876 @@
     if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
       // create all bonds
-      for (int i=0;i<mol->NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
-        OtherFather = mol->ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
+      for (BondList::const_iterator BondRunner = FatherOfRunner->ListOfBonds.begin(); BondRunner != FatherOfRunner->ListOfBonds.end(); (++BondRunner)) {
+        OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
 //        *out << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
         if (SonList[OtherFather->nr] != NULL) {
@@ -887 +884 @@
 //            *out << Verbose(3) << "Adding Bond: ";
 //            *out <<
-            Leaf->AddBond(Runner, SonList[OtherFather->nr], mol->ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree);
+            Leaf->AddBond(Runner, SonList[OtherFather->nr], (*BondRunner)->BondDegree);
 //            *out << "." << endl;
             //NumBonds[Runner->nr]++;
@@ -898 +895 @@
 #ifdef ADDHYDROGEN
           //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
-          if(!Leaf->AddHydrogenReplacementAtom(out, mol->ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, mol->ListOfBondsPerAtom[FatherOfRunner->nr],mol->NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem))
+          if(!Leaf->AddHydrogenReplacementAtom(out, (*BondRunner), Runner, FatherOfRunner, OtherFather, IsAngstroem))
             exit(1);
 #endif
-          //NumBonds[Runner->nr] += ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree;
+          //NumBonds[Runner->nr] += Binder->BondDegree;
         }
       }
@@ -936 +933 @@
   CreateInducedSubgraphOfFragment(out, this, Leaf, SonList, IsAngstroem);
 
-  Leaf->CreateListOfBondsPerAtom(out);
   //Leaflet->Leaf->ScanForPeriodicCorrection(out);
   Free(&SonList);
@@ -981 +977 @@
     if (bit) {  // if bit is set, we add this bond partner
       OtherWalker = BondsSet[j]->rightatom;  // rightatom is always the one more distant, i.e. the one to add
-      //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][CurrentCombination] << ", checking on " << *OtherWalker << "." << endl;
+      //*out << Verbose(1+verbosity) << "Current Bond is " << BondsSet[j] << ", checking on " << *OtherWalker << "." << endl;
       *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
       TestKeySetInsert = FragmentSet->insert(OtherWalker->nr);
@@ -1089 +1085 @@
   NumCombinations = 1 << SetDimension;
 
-  // Hier muessen von 1 bis NumberOfBondsPerAtom[Walker->nr] alle Kombinationen
-  // von Endstuecken (aus den Bonds) hinzugefuegt werden und fuer verbleibende ANOVAOrder
-  // rekursiv GraphCrawler in der naechsten Ebene aufgerufen werden
+  // here for all bonds of Walker all combinations of end pieces (from the bonds)
+  // have to be added and for the remaining ANOVA order GraphCrawler be called
+  // recursively for the next level
 
   *out << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
@@ -1258 +1254 @@
   atom *OtherWalker = NULL;
   atom *Predecessor = NULL;
+  bond *CurrentEdge = NULL;
   bond *Binder = NULL;
-  bond *CurrentEdge = NULL;
   int RootKeyNr = FragmentSearch.Root->GetTrueFather()->nr;
   int RemainingWalkers = -1;
@@ -1285 +1281 @@
       // go through all its bonds
       *out << Verbose(1) << "Going through all bonds of Walker." << endl;
-      for (int i=0;i<mol->NumberOfBondsPerAtom[AtomKeyNr];i++) {
-        Binder = mol->ListOfBondsPerAtom[AtomKeyNr][i];
-        OtherWalker = Binder->GetOtherAtom(Walker);
+      for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+        OtherWalker = (*Runner)->GetOtherAtom(Walker);
         if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
   #ifdef ADDHYDROGEN
@@ -1293 +1288 @@
   #endif
                                                               ) {  // skip hydrogens and restrict to fragment
-          *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
+          *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *(*Runner) << "." << endl;
           // set the label if not set (and push on root stack as well)
           if ((OtherWalker != Predecessor) && (OtherWalker->GetTrueFather()->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
@@ -1712 +1707 @@
         ColorList[Walker->nr] = black;    // mark as explored
         Walker->x.AddVector(&Translationvector); // translate
-        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {  // go through all binding partners
-          if (ListOfBondsPerAtom[Walker->nr][i] != Binder) {
-            OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+          if ((*Runner) != Binder) {
+            OtherWalker = (*Runner)->GetOtherAtom(Walker);
             if (ColorList[OtherWalker->nr] == white) {
               AtomStack->Push(OtherWalker); // push if yet unexplored

src/molecule_geometry.cpp

@@ -258 +258 @@
 {
   atom *Walker = start;
-  bond *Binder = NULL;
   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
   double tmp;
@@ -275 +274 @@
         Testvector.InverseMatrixMultiplication(matrix);
         Translationvector.Zero();
-        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
-          Binder = ListOfBondsPerAtom[Walker->nr][i];
-          if (Walker->nr < Binder->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
+        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+         if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
             for (int j=0;j<NDIM;j++) {
-              tmp = Walker->x.x[j] - Binder->GetOtherAtom(Walker)->x.x[j];
+              tmp = Walker->x.x[j] - (*Runner)->GetOtherAtom(Walker)->x.x[j];
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
-                cout << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *Binder << " has to be shifted due to " << tmp << "." << endl;
+                cout << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl;
                 if (tmp > 0)
                   Translationvector.x[j] -= 1.;
@@ -298 +296 @@
 #ifdef ADDHYDROGEN
         // now also change all hydrogens
-        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
-          Binder = ListOfBondsPerAtom[Walker->nr][i];
-          if (Binder->GetOtherAtom(Walker)->type->Z == 1) {
-            Testvector.CopyVector(&Binder->GetOtherAtom(Walker)->x);
+        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+          if ((*Runner)->GetOtherAtom(Walker)->type->Z == 1) {
+            Testvector.CopyVector(&(*Runner)->GetOtherAtom(Walker)->x);
             Testvector.InverseMatrixMultiplication(matrix);
             Testvector.AddVector(&Translationvector);

src/molecule_graph.cpp

@@ -48 +48 @@
     AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
   }
-
-  CreateListOfBondsPerAtom(out);
 };
 
@@ -64 +62 @@
  *  -# put each atom into its corresponding cell
  *  -# go through every cell, check the atoms therein against all possible bond partners in the 27 adjacent cells, add bond if true
- *  -# create the list of bonds via CreateListOfBondsPerAtom()
  *  -# correct the bond degree iteratively (single->double->triple bond)
  *  -# finally print the bond list to \a *out if desired
@@ -148 +145 @@
     *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << BondDistance << "." << endl;
 
-    // create the adjacency list per atom
-    CreateListOfBondsPerAtom(out);
-
     // correct bond degree by comparing valence and bond degree
     CorrectBondDegree(out);
 
     // output bonds for debugging (if bond chain list was correctly installed)
-    OutputBondsList(out);
+    ActOnAllAtoms( &atom::OutputBondOfAtom, out );
   } else
     *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
@@ -173 +167 @@
   }
   *out << endl;
-};
-
-/** Sums up the number of bonds times bond::BondDegree the atom with index \a nr has.
- * \param nr indexof atom
- * \return number of bonds each weighted with its bond::BondDegree
- */
-int molecule::CountAtomsBonds(int nr)
-{
-  int NoBonds = 0;
-  for(int j=0;j<NumberOfBondsPerAtom[nr];j++)
-    NoBonds += ListOfBondsPerAtom[nr][j]->BondDegree;
-  return NoBonds;
 };
 
@@ -198 +180 @@
 {
   int No = 0;
-  int NoBonds = 0;
-  int CandidateBondNo = 0;
-  int FalseBondDegree = 0;
-  atom *Walker = NULL;
-  atom *Candidate = NULL;
-  atom *OtherWalker = NULL;
 
   if (BondCount != 0) {
-    NoCyclicBonds = 0;
-    *out << Verbose(1) << "Correcting Bond degree of each bond ... ";
+    *out << Verbose(1) << "Correcting Bond degree of each bond ... " << endl;
     do {
-      No = 0; // No acts as breakup flag (if 1 we still continue)
-      Walker = start;
-      while (Walker->next != end) { // go through every atom
-        Walker = Walker->next;
-
-        NoBonds = CountAtomsBonds(Walker->nr);
-        *out << Verbose(3) << "Walker " << *Walker << ": " << (int)Walker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-        if ((int)(Walker->type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
-          Candidate = NULL;
-          CandidateBondNo = -1;
-          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through each of its bond partners
-            OtherWalker = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-            NoBonds = CountAtomsBonds(OtherWalker->nr);
-            *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
-            if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
-              if ((Candidate == NULL) || (NumberOfBondsPerAtom[Candidate->nr] > NumberOfBondsPerAtom[OtherWalker->nr])) { // pick the one with fewer number of bonds first
-                Candidate = OtherWalker;
-                CandidateBondNo = i;
-                *out << Verbose(3) << "New candidate is " << *Candidate << "." << endl;
-              }
-            }
-          }
-          if ((Candidate != NULL) && (CandidateBondNo != -1)) {
-            ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
-            *out << Verbose(2) << "Increased bond degree for bond " << *ListOfBondsPerAtom[Walker->nr][CandidateBondNo] << "." << endl;
-          } else {
-            *out << Verbose(2) << "Could not find correct degree for atom " << *Walker << "." << endl;
-            FalseBondDegree++;
-          }
-        }
-      }
+      No = SumPerAtom( &atom::CorrectBondDegree, out );
     } while (No);
     *out << " done." << endl;
@@ -245 +190 @@
     *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
   }
-  *out << FalseBondDegree << " bonds could not be corrected." << endl;
-
-  return (FalseBondDegree);
+  *out << No << " bonds could not be corrected." << endl;
+
+  return (No);
 };
 
@@ -257 +202 @@
 int molecule::CountCyclicBonds(ofstream *out)
 {
-  int No = 0;
+  NoCyclicBonds = 0;
   int *MinimumRingSize = NULL;
   MoleculeLeafClass *Subgraphs = NULL;
@@ -275 +220 @@
     Binder = Binder->next;
     if (Binder->Cyclic)
-      No++;
+      NoCyclicBonds++;
   }
   delete(BackEdgeStack);
-  return No;
+  return NoCyclicBonds;
 };
 
@@ -459 +404 @@
   }
   // set cyclic bond criterium on "same LP" basis
-  Binder = first;
+  CyclicBondAnalysis();
+
+  OutputGraphInfoPerAtom(out);
+
+  OutputGraphInfoPerBond(out);
+
+  // free all and exit
+  delete(AtomStack);
+  *out << Verbose(0) << "End of DepthFirstSearchAnalysis" << endl;
+  return SubGraphs;
+};
+
+/** Scans through all bonds and set bond::Cyclic to true where atom::LowpointNr of both ends is equal: LP criterion.
+ */
+void molecule::CyclicBondAnalysis()
+{
+  NoCyclicBonds = 0;
+  bond *Binder = first;
   while(Binder->next != last) {
     Binder = Binder->next;
@@ -467 +429 @@
     }
   }
-
-
+};
+
+/** Output graph information per atom.
+ * \param *out output stream
+ */
+void molecule::OutputGraphInfoPerAtom(ofstream *out)
+{
   *out << Verbose(1) << "Final graph info for each atom is:" << endl;
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    *out << Verbose(2) << "Atom " << Walker->Name << " is " << ((Walker->SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
-    OutputComponentNumber(out, Walker);
-    *out << " with Lowpoint " << Walker->LowpointNr << " and Graph Nr. " << Walker->GraphNr << "." << endl;
-  }
-
+  ActOnAllAtoms( &atom::OutputGraphInfo, out );
+};
+
+/** Output graph information per bond.
+ * \param *out output stream
+ */
+void molecule::OutputGraphInfoPerBond(ofstream *out)
+{
   *out << Verbose(1) << "Final graph info for each bond is:" << endl;
-  Binder = first;
+  bond *Binder = first;
   while(Binder->next != last) {
     Binder = Binder->next;
     *out << Verbose(2) << ((Binder->Type == TreeEdge) ? "TreeEdge " : "BackEdge ") << *Binder << ": <";
     *out << ((Binder->leftatom->SeparationVertex) ? "SP," : "") << "L" << Binder->leftatom->LowpointNr << " G" << Binder->leftatom->GraphNr << " Comp.";
-    OutputComponentNumber(out, Binder->leftatom);
+    Binder->leftatom->OutputComponentNumber(out);
     *out << " ===  ";
     *out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
-    OutputComponentNumber(out, Binder->rightatom);
+    Binder->rightatom->OutputComponentNumber(out);
     *out << ">." << endl;
     if (Binder->Cyclic) // cyclic ??
       *out << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl;
   }
-
-  // free all and exit
-  delete(AtomStack);
-  *out << Verbose(0) << "End of DepthFirstSearchAnalysis" << endl;
-  return SubGraphs;
 };
 
@@ -547 +509 @@
       Walker = BFSStack->PopFirst();
       *out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
-      for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-        Binder = ListOfBondsPerAtom[Walker->nr][i];
-        if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
-          OtherAtom = Binder->GetOtherAtom(Walker);
+      for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+        if ((*Runner) != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
+          OtherAtom = (*Runner)->GetOtherAtom(Walker);
 #ifdef ADDHYDROGEN
           if (OtherAtom->type->Z != 1) {
@@ -662 +623 @@
           Walker = BFSStack->PopFirst();
           //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
-          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-            Binder = ListOfBondsPerAtom[Walker->nr][i];
-            if ((Binder != BackEdge) || (NumberOfBondsPerAtom[Walker->nr] == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check
-              OtherAtom = Binder->GetOtherAtom(Walker);
+          for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+            if (((*Runner) != BackEdge) || (Walker->ListOfBonds.size() == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check
+              OtherAtom = (*Runner)->GetOtherAtom(Walker);
               //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
               if (ColorList[OtherAtom->nr] == white) {
@@ -717 +677 @@
 void molecule::SetNextComponentNumber(atom *vertex, int nr)
 {
-  int i=0;
+  size_t i=0;
   if (vertex != NULL) {
-    for(;i<NumberOfBondsPerAtom[vertex->nr];i++) {
+    for(;i<vertex->ListOfBonds.size();i++) {
       if (vertex->ComponentNr[i] == -1) {   // check if not yet used
         vertex->ComponentNr[i] = nr;
@@ -727 +687 @@
         break;  // breaking here will not cause error!
     }
-    if (i == NumberOfBondsPerAtom[vertex->nr])
+    if (i == vertex->ListOfBonds.size())
       cerr << "Error: All Component entries are already occupied!" << endl;
   } else
       cerr << "Error: Given vertex is NULL!" << endl;
-};
-
-/** Output a list of flags, stating whether the bond was visited or not.
- * \param *out output stream for debugging
- */
-void molecule::OutputComponentNumber(ofstream *out, atom *vertex)
-{
-  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
-    *out << vertex->ComponentNr[i] << "  ";
 };
 
@@ -751 +702 @@
     if (Walker->ComponentNr != NULL)
       Free(&Walker->ComponentNr);
-    Walker->ComponentNr = Malloc<int>(NumberOfBondsPerAtom[Walker->nr], "molecule::InitComponentNumbers: *Walker->ComponentNr");
-    for (int i=NumberOfBondsPerAtom[Walker->nr];i--;)
+    Walker->ComponentNr = Malloc<int>(Walker->ListOfBonds.size()+1, "molecule::InitComponentNumbers: *Walker->ComponentNr");
+    for (int i=Walker->ListOfBonds.size()+1;i--;)
       Walker->ComponentNr[i] = -1;
   }
@@ -763 +714 @@
 bond * molecule::FindNextUnused(atom *vertex)
 {
-  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
-    if (ListOfBondsPerAtom[vertex->nr][i]->IsUsed() == white)
-      return(ListOfBondsPerAtom[vertex->nr][i]);
+  for (BondList::const_iterator Runner = vertex->ListOfBonds.begin(); Runner != vertex->ListOfBonds.end(); (++Runner))
+    if ((*Runner)->IsUsed() == white)
+      return((*Runner));
   return NULL;
 };
@@ -813 +764 @@
 {
   ofstream AdjacencyFile;
-  atom *Walker = NULL;
   stringstream line;
   bool status = true;
@@ -821 +771 @@
   *out << Verbose(1) << "Saving adjacency list ... ";
   if (AdjacencyFile != NULL) {
-    Walker = start;
-    while(Walker->next != end) {
-      Walker = Walker->next;
-      AdjacencyFile << Walker->nr << "\t";
-      for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
-        AdjacencyFile << ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker)->nr << "\t";
-      AdjacencyFile << endl;
-    }
+    ActOnAllAtoms( &atom::OutputAdjacency, &AdjacencyFile );
     AdjacencyFile.close();
     *out << Verbose(1) << "done." << endl;
@@ -850 +793 @@
   stringstream filename;
   bool status = true;
+  atom *Walker = NULL;
   char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
 
@@ -859 +803 @@
     int NonMatchNumber = 0;   // will number of atoms with differing bond structure
     int *CurrentBonds = Malloc<int>(8, "molecule::CheckAdjacencyFileAgainstMolecule - CurrentBonds"); // contains parsed bonds of current atom
-    int CurrentBondsOfAtom;
+    size_t CurrentBondsOfAtom;
 
     // Parse the file line by line and count the bonds
@@ -871 +815 @@
       // parse into structure
       if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
+        Walker = ListOfAtoms[AtomNr];
         while (!line.eof())
           line >> CurrentBonds[ ++CurrentBondsOfAtom ];
         // compare against present bonds
         //cout << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
-        if (CurrentBondsOfAtom == NumberOfBondsPerAtom[AtomNr]) {
-          for(int i=0;i<NumberOfBondsPerAtom[AtomNr];i++) {
-            int id = ListOfBondsPerAtom[AtomNr][i]->GetOtherAtom(ListOfAtoms[AtomNr])->nr;
-            int j = 0;
+        if (CurrentBondsOfAtom == Walker->ListOfBonds.size()) {
+          for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+            int id = (*Runner)->GetOtherAtom(Walker)->nr;
+            size_t j = 0;
             for (;(j<CurrentBondsOfAtom) && (CurrentBonds[j++] != id);); // check against all parsed bonds
             if (CurrentBonds[j-1] != id) { // no match ? Then mark in ListOfAtoms
@@ -891 +836 @@
           //out << endl;
         } else {
-          *out << "Number of bonds for Atom " << *ListOfAtoms[AtomNr] << " does not match, parsed " << CurrentBondsOfAtom << " against " << NumberOfBondsPerAtom[AtomNr] << "." << endl;
+          *out << "Number of bonds for Atom " << *Walker << " does not match, parsed " << CurrentBondsOfAtom << " against " << Walker->ListOfBonds.size() << "." << endl;
           status = false;
         }
@@ -937 +882 @@
     Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
     if (Walker != NULL) // if this Walker exists in the subgraph ...
-      for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
-        OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-        if (OtherAtom == ListOfLocalAtoms[Binder->rightatom->nr]) { // found the bond
-          LocalStack->Push(ListOfBondsPerAtom[Walker->nr][i]);
-          *out << Verbose(3) << "Found local edge " << *(ListOfBondsPerAtom[Walker->nr][i]) << "." << endl;
+      for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+        OtherAtom = (*Runner)->GetOtherAtom(Walker);
+        if (OtherAtom == ListOfLocalAtoms[(*Runner)->rightatom->nr]) { // found the bond
+          LocalStack->Push((*Runner));
+          *out << Verbose(3) << "Found local edge " << *(*Runner) << "." << endl;
           break;
         }
@@ -973 +918 @@
   class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
   atom *Walker = NULL, *OtherAtom = NULL;
-  bond *Binder = NULL;
 
   // add Root if not done yet
@@ -999 +943 @@
     // followed by n+1 till top of stack.
     Walker = AtomStack->PopFirst(); // pop oldest added
-    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl;
-    for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-      Binder = ListOfBondsPerAtom[Walker->nr][i];
-      if (Binder != NULL) { // don't look at bond equal NULL
-        OtherAtom = Binder->GetOtherAtom(Walker);
-        *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
+    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << Walker->ListOfBonds.size() << " bonds." << endl;
+    for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+      if ((*Runner) != NULL) { // don't look at bond equal NULL
+        OtherAtom = (*Runner)->GetOtherAtom(Walker);
+        *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *(*Runner) << "." << endl;
         if (ColorList[OtherAtom->nr] == white) {
-          if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
+          if ((*Runner) != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
             ColorList[OtherAtom->nr] = lightgray;
           PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
           ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
           *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
-          if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond))) ) { // Check for maximum distance
+          if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && ((*Runner) != Bond))) ) { // Check for maximum distance
             *out << Verbose(3);
             if (AddedAtomList[OtherAtom->nr] == NULL) { // add if it's not been so far
               AddedAtomList[OtherAtom->nr] = Mol->AddCopyAtom(OtherAtom);
               *out << "Added OtherAtom " << OtherAtom->Name;
-              AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-              AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-              AddedBondList[Binder->nr]->Type = Binder->Type;
-              *out << " and bond " << *(AddedBondList[Binder->nr]) << ", ";
+              AddedBondList[(*Runner)->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], (*Runner));
+              *out << " and bond " << *(AddedBondList[(*Runner)->nr]) << ", ";
             } else {  // this code should actually never come into play (all white atoms are not yet present in BondMolecule, that's why they are white in the first place)
               *out << "Not adding OtherAtom " << OtherAtom->Name;
-              if (AddedBondList[Binder->nr] == NULL) {
-                AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-                AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-                AddedBondList[Binder->nr]->Type = Binder->Type;
-                *out << ", added Bond " << *(AddedBondList[Binder->nr]);
+              if (AddedBondList[(*Runner)->nr] == NULL) {
+                AddedBondList[(*Runner)->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], (*Runner));
+                *out << ", added Bond " << *(AddedBondList[(*Runner)->nr]);
               } else
                 *out << ", not added Bond ";
@@ -1033 +972 @@
             AtomStack->Push(OtherAtom);
           } else { // out of bond order, then replace
-            if ((AddedAtomList[OtherAtom->nr] == NULL) && (Binder->Cyclic))
+            if ((AddedAtomList[OtherAtom->nr] == NULL) && ((*Runner)->Cyclic))
               ColorList[OtherAtom->nr] = white; // unmark if it has not been queued/added, to make it available via its other bonds (cyclic)
-            if (Binder == Bond)
+            if ((*Runner) == Bond)
               *out << Verbose(3) << "Not Queueing, is the Root bond";
             else if (ShortestPathList[OtherAtom->nr] >= BondOrder)
               *out << Verbose(3) << "Not Queueing, is out of Bond Count of " << BondOrder;
-            if (!Binder->Cyclic)
+            if (!(*Runner)->Cyclic)
               *out << ", is not part of a cyclic bond, saturating bond with Hydrogen." << endl;
-            if (AddedBondList[Binder->nr] == NULL) {
+            if (AddedBondList[(*Runner)->nr] == NULL) {
               if ((AddedAtomList[OtherAtom->nr] != NULL)) { // .. whether we add or saturate
-                AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-                AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-                AddedBondList[Binder->nr]->Type = Binder->Type;
+                AddedBondList[(*Runner)->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], (*Runner));
               } else {
 #ifdef ADDHYDROGEN
-                if (!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem))
+                if (!Mol->AddHydrogenReplacementAtom(out, (*Runner), AddedAtomList[Walker->nr], Walker, OtherAtom, IsAngstroem))
                   exit(1);
 #endif
@@ -1057 +994 @@
           *out << Verbose(3) << "Not Adding, has already been visited." << endl;
           // This has to be a cyclic bond, check whether it's present ...
-          if (AddedBondList[Binder->nr] == NULL) {
-            if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
-              AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
-              AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
-              AddedBondList[Binder->nr]->Type = Binder->Type;
+          if (AddedBondList[(*Runner)->nr] == NULL) {
+            if (((*Runner) != Bond) && ((*Runner)->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
+              AddedBondList[(*Runner)->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], (*Runner));
             } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
 #ifdef ADDHYDROGEN
-              if(!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem))
+              if(!Mol->AddHydrogenReplacementAtom(out, (*Runner), AddedAtomList[Walker->nr], Walker, OtherAtom, IsAngstroem))
                 exit(1);
 #endif
@@ -1080 +1015 @@
   delete(AtomStack);
 };
+
+/** Adds a bond as a copy to a given one
+ * \param *left leftatom of new bond
+ * \param *right rightatom of new bond
+ * \param *CopyBond rest of fields in bond are copied from this
+ * \return pointer to new bond
+ */
+bond * molecule::CopyBond(atom *left, atom *right, bond *CopyBond)
+{
+  bond *Binder = AddBond(left, right, CopyBond->BondDegree);
+  Binder->Cyclic = CopyBond->Cyclic;
+  Binder->Type = CopyBond->Type;
+  return Binder;
+};
+
 
 /** Adds bond structure to this molecule from \a Father molecule.
@@ -1122 +1072 @@
         status = false;
       } else {
-        for (int i=0;i<Father->NumberOfBondsPerAtom[Walker->nr];i++) {
-          OtherAtom = Father->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+          OtherAtom = (*Runner)->GetOtherAtom(Walker);
           if (ParentList[OtherAtom->nr] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
-            *out << Verbose(4) << "Endpoints of Bond " << Father->ListOfBondsPerAtom[Walker->nr][i] << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
-            AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], Father->ListOfBondsPerAtom[Walker->nr][i]->BondDegree);
+            *out << Verbose(4) << "Endpoints of Bond " << (*Runner) << " are both present: " << ParentList[Walker->nr]->Name << " and " << ParentList[OtherAtom->nr]->Name << "." << endl;
+            AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], (*Runner)->BondDegree);
           }
         }
@@ -1164 +1114 @@
       for(KeySet::iterator runners = Fragment->begin(); runners != Fragment->end(); runners++) {
         Walker2 = FindAtom(*runners);
-        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
-          if (ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker) == Walker2) {
+        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
+          if ((*Runner)->GetOtherAtom(Walker) == Walker2) {
             BondStatus = true;
             break;

src/molecule_template.hpp

@@ -153 +153 @@
 };
 
+// ========================= Summing over each Atoms =================================== //
+
+// zero arguments
+template <typename res> res molecule::SumPerAtom(res (atom::*f)() )
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)();
+  }
+  return result;
+};
+template <typename res> res molecule::SumPerAtom(res (atom::*f)() const )
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)();
+  }
+  return result;
+};
+template <typename res> res molecule::SumPerAtom(res (atom::*f)() ) const
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)();
+  }
+  return result;
+};
+template <typename res> res molecule::SumPerAtom(res (atom::*f)() const ) const
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)();
+  }
+  return result;
+};
+// one argument
+template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T), T t )
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)(t);
+  }
+  return result;
+};
+template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T) const, T t )
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)(t);
+  }
+  return result;
+};
+template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T), T t ) const
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)(t);
+  }
+  return result;
+};
+template <typename T, typename res> res molecule::SumPerAtom(res (atom::*f)(T) const, T t ) const
+{
+  res result = 0;
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    result += (Walker->*f)(t);
+  }
+  return result;
+};
+
+
 // ================== Acting with each Atoms on same molecule ========================== //
 

src/moleculelist.cpp

@@ -376 +376 @@
   atom *Walker = NULL;
   atom *Runner = NULL;
+  bond *Binder = NULL;
   double ***FitConstant = NULL, **correction = NULL;
   int a, b;
@@ -482 +483 @@
     while (Walker->next != (*ListRunner)->end) {
       Walker = Walker->next;
-      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
+      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
       if ((Walker->type->Z == 1) && ((Walker->father == NULL)
           || (Walker->father->type->Z != 1))) { // if it's a hydrogen
@@ -488 +489 @@
         while (Runner->next != (*ListRunner)->end) {
           Runner = Runner->next;
-          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
+          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
           // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
-          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && ((*ListRunner)->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != (*ListRunner)->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
+          Binder = *(Runner->ListOfBonds.begin());
+          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (Binder->GetOtherAtom(Runner) != Binder->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
             // 4. evaluate the morse potential for each matrix component and add up
             distance = Runner->x.Distance(&Walker->x);
@@ -812 +814 @@
 bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
 {
-  atom *Walker = NULL, *OtherWalker = NULL;
-  bond *Binder = NULL;
+  atom *Walker = NULL;
+  atom *OtherWalker = NULL;
+  atom *Father = NULL;
   bool status = true;
   int AtomNo;
@@ -830 +833 @@
     while (Walker->next != Leaf->end) {
       Walker = Walker->next;
-      AtomNo = Walker->GetTrueFather()->nr; // global id of the current walker
-      for (int i = 0; i < reference->NumberOfBondsPerAtom[AtomNo]; i++) { // go through father's bonds and copy them all
-        Binder = reference->ListOfBondsPerAtom[AtomNo][i];
-        OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
+      Father = Walker->GetTrueFather();
+      AtomNo = Father->nr; // global id of the current walker
+      for (BondList::const_iterator Runner = Father->ListOfBonds.begin(); Runner != Father->ListOfBonds.end(); (++Runner)) {
+        OtherWalker = ListOfLocalAtoms[FragmentCounter][(*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
         if (OtherWalker != NULL) {
           if (OtherWalker->nr > Walker->nr)
-            Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
+            Leaf->AddBond(Walker, OtherWalker, (*Runner)->BondDegree);
         } else {
-          *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
+          *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << (*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
           status = false;
         }
       }
     }
-    Leaf->CreateListOfBondsPerAtom(out);
     FragmentCounter++;
     if (next != NULL)

src/unittests/Makefile.am

@@ -1 +1 @@
 INCLUDES = -I$(top_srcdir)/src
 
-noinst_PROGRAMS =  ActOnAllUnitTest AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest MemoryAllocatorUnitTest MemoryUsageObserverUnitTest VectorUnitTest
+noinst_PROGRAMS =  ActOnAllUnitTest AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest ListOfBondsUnitTest MemoryAllocatorUnitTest MemoryUsageObserverUnitTest VectorUnitTest
 
-TESTS = ActOnAllUnitTest AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest VectorUnitTest
+TESTS = ActOnAllUnitTest AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest ListOfBondsUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest VectorUnitTest
 check_PROGRAMS = $(TESTS)
 
@@ -26 +26 @@
 AnalysisPairCorrelationUnitTest_LDADD = ../libmolecuilder.a
 
+ListOfBondsUnitTest_SOURCES = listofbondsunittest.cpp listofbondsunittest.hpp
+ListOfBondsUnitTest_CXXFLAGS = $(CPPUNIT_CFLAGS)
+ListOfBondsUnitTest_LDFLAGS = $(CPPUNIT_LIBS) -ldl
+ListOfBondsUnitTest_LDADD = ../libmolecuilder.a
+
 VectorUnitTest_SOURCES = defs.hpp helpers.hpp leastsquaremin.hpp memoryallocator.hpp memoryusageobserver.hpp vectorunittest.cpp vectorunittest.hpp vector.hpp verbose.hpp 
 VectorUnitTest_CXXFLAGS = $(CPPUNIT_CFLAGS)