Changeset 3738f0 for src/molecule_graph.cpp

Timestamp:

Mar 1, 2011, 12:49:28 PM (14 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:

0cbad2

Parents:

ce5f05

git-author:

Frederik Heber <heber@…> (02/26/11 00:10:05)

git-committer:

Frederik Heber <heber@…> (03/01/11 12:49:28)

Message:

Moved molecule::CreateAdjacencyList over to class BondGraph.

to make this possible we had to:

• enhance AtomSet such that individual const_ and iterator are templated.
• FIX: FragmentationCreateAdjacencyAction needs to call BondGraph::SetMaxDistanceToMaxOfCovalentRadii() before getMaxDistance() is valid.
• new function molecule::getAtomSet() which copies the internal atomSet to a vector AtomSetMixin and delivers. This is needed as PointCloudAdaptor needs a reference to a AtomSetMixin and we don't want to expose the internal set.

other changes:

• the following Actions do not create adjacency beforehand anymore:
  • DepthFirstSearchAction
  • MoleculeSaveAdjacencyAction
  • MoleculeSaveBondsAction
• changes to BondGraph:
  • removed CreateAdjacencyForMolecule()
  • removed SetMaxDistanceToMaxOfCovalentRadii()
  • new function CreateAdjacency on LinkedCell reference
• new BondGraph functions that work on AtomSetMixin reference: CreateAdjacency(), cleanAdjacencyList(), CorrectBondDegree(), resetBondDegree(), calculateBondDegree()

TESTFIXES:

• the regression test for all Actions mentioned above that don't create adjacency themselves anymore needed to be prepended with --select-all-atoms --create-adjacency.

File:

• src/molecule_graph.cpp (modified) (2 diffs)

src/molecule_graph.cpp

@@ -111 +111 @@
 }
 
-/** Creates an adjacency list of the molecule.
- * Generally, we use the CSD approach to bond recognition, that is the the distance
- * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
- * a threshold t = 0.4 Angstroem.
- * To make it O(N log N) the function uses the linked-cell technique as follows:
- * The procedure is step-wise:
- *  -# Remove every bond in list
- *  -# Count the atoms in the molecule with CountAtoms()
- *  -# partition cell into smaller linked cells of size \a bonddistance
- *  -# 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
- *  -# correct the bond degree iteratively (single->double->triple bond)
- *  -# finally print the bond list to \a *out if desired
- * \param bonddistance length of linked cells (i.e. maximum minimal length checked)
- * \param IsAngstroem whether coordinate system is gauged to Angstroem or Bohr radii
- * \param *minmaxdistance function to give upper and lower bound on whether particle is bonded to some other
- * \param *BG BondGraph with the member function above or NULL, if just standard covalent should be used.
- */
-void molecule::CreateAdjacencyList(
-    double bonddistance,
-    bool IsAngstroem,
-    void (BondGraph::*minmaxdistance)(const BondedParticle * const , const BondedParticle * const , double &, double &, bool),
-    BondGraph *BG)
-{
-  atom *Walker = NULL;
-  atom *OtherWalker = NULL;
-  int n[NDIM];
-  double MinDistance, MaxDistance;
-  LinkedCell *LC = NULL;
-  bool free_BG = false;
-  Box &domain = World::getInstance().getDomain();
-
-  DoLog(0) && (Log() << Verbose(0) << "Begin of CreateAdjacencyList." << endl);
-  // remove every bond from the list
-  for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
-    BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
-    for(BondList::iterator BondRunner = ListOfBonds.begin();
-        !ListOfBonds.empty();
-        BondRunner = ListOfBonds.begin())
-      if ((*BondRunner)->leftatom == *AtomRunner)
-        delete((*BondRunner));
-  }
-
-  // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
-  DoLog(1) && (Log() << Verbose(1) << "AtomCount " << getAtomCount() << " and bonddistance is " << bonddistance << "." << endl);
-
-  if ((getAtomCount() > 1) && (bonddistance > 0.1)) {
-    DoLog(2) && (Log() << Verbose(2) << "Creating Linked Cell structure ... " << endl);
-    PointCloudAdaptor<molecule> cloud(this, name);
-    LC = new LinkedCell(cloud, bonddistance);
-
-    // create a list to map Tesselpoint::Nr to atom *
-    DoLog(2) && (Log() << Verbose(2) << "Creating TesselPoint to atom map ... " << endl);
-
-    // set numbers for atoms that can later be used
-    int i=0;
-    for(internal_iterator iter = atoms.begin();iter!= atoms.end(); ++iter){
-      (*iter)->setNr(i++);
-    }
-
-    // 3a. go through every cell
-    DoLog(2) && (Log() << Verbose(2) << "Celling ... " << endl);
-    for (LC->n[0] = 0; LC->n[0] < LC->N[0]; LC->n[0]++)
-      for (LC->n[1] = 0; LC->n[1] < LC->N[1]; LC->n[1]++)
-        for (LC->n[2] = 0; LC->n[2] < LC->N[2]; LC->n[2]++) {
-          const TesselPointSTLList *List = LC->GetCurrentCell();
-          Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
-          if (List != NULL) {
-            for (TesselPointSTLList::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-              Walker = dynamic_cast<atom*>(*Runner);
-              ASSERT(Walker,"Tesselpoint that was not an atom retrieved from LinkedNode");
-              Log() << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
-              // 3c. check for possible bond between each atom in this and every one in the 27 cells
-              for (n[0] = -1; n[0] <= 1; n[0]++)
-                for (n[1] = -1; n[1] <= 1; n[1]++)
-                  for (n[2] = -1; n[2] <= 1; n[2]++) {
-                    const TesselPointSTLList *OtherList = LC->GetRelativeToCurrentCell(n);
-                    if (OtherList != NULL) {
-                      Log() << Verbose(2) << "Current relative cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
-                      for (TesselPointSTLList::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
-                        if ((*OtherRunner)->getNr() > Walker->getNr()) {
-                          OtherWalker = dynamic_cast<atom*>(*OtherRunner);
-                          ASSERT(OtherWalker,"TesselPoint that was not an atom retrieved from LinkedNode");
-                          (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
-                          const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
-                          Log() << Verbose(1) << "Checking distance " << distance << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
-                          const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
-                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
-                          if (OtherWalker->father->getNr() > Walker->father->getNr()) {
-                            if (status) { // create bond if distance is smaller
-                              Log() << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
-                              AddBond(Walker->father, OtherWalker->father, 1); // also increases molecule::BondCount
-                            } else {
-                              Log() << Verbose(1) << "Not Adding: distance too great." << endl;
-                            }
-                          } else {
-                            Log() << Verbose(1) << "Not Adding: Wrong order of labels." << endl;
-                          }
-                        }
-                      }
-                    }
-                  }
-            }
-          }
-        }
-    delete (LC);
-    DoLog(1) && (Log() << Verbose(1) << "I detected " << getBondCount() << " bonds in the molecule." << endl);
-
-    // correct bond degree by comparing valence and bond degree
-    DoLog(2) && (Log() << Verbose(2) << "Correcting bond degree ... " << endl);
-    CorrectBondDegree();
-
-    // output bonds for debugging (if bond chain list was correctly installed)
-    for(molecule::internal_iterator iter = atoms.begin(); iter != atoms.end(); ++iter)
-      (*iter)->OutputBondOfAtom((std::ostream &)std::cout);
-  } else
-    DoLog(1) && (Log() << Verbose(1) << "AtomCount is " << getAtomCount() << ", thus no bonds, no connections!." << endl);
-  DoLog(0) && (Log() << Verbose(0) << "End of CreateAdjacencyList." << endl);
-  if (free_BG)
-    delete(BG);
-}
-;
-
 /** Checks for presence of bonds within atom list.
  * TODO: more sophisticated check for bond structure (e.g. connected subgraph, ...)
@@ -265 +142 @@
 }
 ;
-
-/** correct bond degree by comparing valence and bond degree.
- * correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
- * iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
- * preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
- * double bonds as was expected.
- * \return number of bonds that could not be corrected
- */
-int molecule::CorrectBondDegree() const
-{
-  int No = 0, OldNo = -1;
-
-  if (getBondCount() != 0) {
-    DoLog(1) && (Log() << Verbose(1) << "Correcting Bond degree of each bond ... " << endl);
-    do {
-      OldNo = No;
-      No=0;
-      BOOST_FOREACH(atom *atom,atoms){
-        No+=atom->CorrectBondDegree();
-      }
-    } while (OldNo != No);
-    DoLog(0) && (Log() << Verbose(0) << " done." << endl);
-  } else {
-    DoLog(1) && (Log() << Verbose(1) << "BondCount is " << getBondCount() << ", no bonds between any of the " << getAtomCount() << " atoms." << endl);
-  }
-  DoLog(0) && (Log() << Verbose(0) << No << " bonds could not be corrected." << endl);
-
-  return (No);
-}