Changeset 826e8c for src


Ignore:
Timestamp:
Feb 22, 2012, 11:29:04 AM (13 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:
153985
Parents:
cbcda6
git-author:
Frederik Heber <heber@…> (02/06/12 17:46:56)
git-committer:
Frederik Heber <heber@…> (02/22/12 11:29:04)
Message:

Introduced new LinkedCell into BondGraph.

  • new convenience functions BondGraph::getDomain(), ::getLinkedCell(), ::getTime() to have World dependency not in header file (circular include).
  • new function BondGraph::getMaxPossibleBondDistance() where one element is given.
  • dropped BondGraph::CreateAdjacency(LC_deprecated) as functionality is better placed directly into templated version directly.
  • rewrote templated BondGraph::CreateAdjacency that now makes use of LinkedCell_View::getAllNeighbours().
  • FIX: BondGraph::CreateAdjacency() created bonds for atom's fathers, not for atoms.
  • TESTFIX: New LinkedCell construct in BondGraph::CreateAdjacency changes order of ids in adjacency file:
    • ids are now in order (unlike to before), hence we simply use the new file in the regression test.
    • this is actually a fault of the test, not of the code as order of ids is not guaranteed in the file anyway.
  • TESTFIX: Removed XFAILs from regression test Graph/SubgraphDissection- BoundaryCondition.
  • TESTFIX: CountBondsUnitTest use two molecules that are on top of each other: This causes mess with changed LinkedCell BondGraph recognition, and we even actually need just one molecule for the testing.
Location:
src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • src/Analysis/unittests/CountBondsUnitTest.cpp

    rcbcda6 r826e8c  
    3333#include "Atom/atom.hpp"
    3434#include "Bond/bond.hpp"
     35#include "CodePatterns/Log.hpp"
    3536#include "Element/element.hpp"
    3637#include "Element/periodentafel.hpp"
    … …  
    6263  CPPUNIT_ASSERT(oxygen != NULL && "could not find element oxygen");
    6364
     65  //setVerbosity(3);
     66
    6467  // construct molecule (water molecule)
    6568  TestMolecule1 = World::getInstance().createMolecule();
    66   CPPUNIT_ASSERT(TestMolecule1 != NULL && "could not create first molecule");
     69  CPPUNIT_ASSERT(TestMolecule1 != NULL && "could not create second molecule");
    6770  Walker = World::getInstance().createAtom();
    6871  CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
    … …  
    8184  TestMolecule1->AddAtom(Walker);
    8285
    83   TestMolecule2 = World::getInstance().createMolecule();
    84   CPPUNIT_ASSERT(TestMolecule2 != NULL && "could not create second molecule");
    85   Walker = World::getInstance().createAtom();
    86   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
    87   Walker->setType(hydrogen);
    88   Walker->setPosition(Vector(-0.2418, 0.9350, 0. ));
    89   TestMolecule2->AddAtom(Walker);
    90   Walker = World::getInstance().createAtom();
    91   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
    92   Walker->setType(hydrogen);
    93   Walker->setPosition(Vector(0.9658, 0., 0. ));
    94   TestMolecule2->AddAtom(Walker);
    95   Walker = World::getInstance().createAtom();
    96   CPPUNIT_ASSERT(Walker != NULL && "could not create atom");
    97   Walker->setType(oxygen);
    98   Walker->setPosition(Vector(0., 0., 0. ));
    99   TestMolecule2->AddAtom(Walker);
    100 
    10186  molecules = World::getInstance().getMolecules();
    10287  CPPUNIT_ASSERT(molecules != NULL && "could not obtain list of molecules");
    10388  molecules->insert(TestMolecule1);
    104   molecules->insert(TestMolecule2);
    105 
    106   // check that TestMolecule was correctly constructed
     89
     90  // check that TestMolecule1 was correctly constructed
    10791  CPPUNIT_ASSERT_EQUAL( TestMolecule1->getAtomCount(), 3 );
    108   CPPUNIT_ASSERT_EQUAL( TestMolecule2->getAtomCount(), 3 );
    10992
    11093  // create a small file with table
    … …  
    11699  BG->CreateAdjacency(Set1);
    117100  CPPUNIT_ASSERT( TestMolecule1->getBondCount() != 0);
    118   World::AtomComposite Set2 = TestMolecule2->getAtomSet();
    119   BG->CreateAdjacency(Set2);
    120   CPPUNIT_ASSERT( TestMolecule2->getBondCount() != 0);
    121101//  TestMolecule1->Output((ofstream *)&cout);
    122102//  TestMolecule1->OutputBondsList(std::cout);
    … …  
    136116void CountBondsTest::BondsOfTwoTest()
    137117{
    138   CPPUNIT_ASSERT_EQUAL( 4 , CountBondsOfTwo(molecules, hydrogen, oxygen) );
     118  CPPUNIT_ASSERT_EQUAL( 2 , CountBondsOfTwo(molecules, hydrogen, oxygen) );
    139119  CPPUNIT_ASSERT_EQUAL( 0 , CountBondsOfTwo(molecules, hydrogen, hydrogen) );
    140120  CPPUNIT_ASSERT_EQUAL( 0 , CountBondsOfTwo(molecules, oxygen, oxygen) );
    … …  
    145125void CountBondsTest::BondsOfThreeTest()
    146126{
    147   CPPUNIT_ASSERT_EQUAL( 2 , CountBondsOfThree(molecules, hydrogen, oxygen, hydrogen) );
     127  CPPUNIT_ASSERT_EQUAL( 1 , CountBondsOfThree(molecules, hydrogen, oxygen, hydrogen) );
    148128  CPPUNIT_ASSERT_EQUAL( 0 , CountBondsOfThree(molecules, oxygen, hydrogen, oxygen) );
    149129};
    150130
    151 void OutputTestMolecule(molecule *mol, const char *name)
     131void OutputTestMolecule1(molecule *mol, const char *name)
    152132{
    153133  ofstream output(name);
    … …  
    166146  Vector Translator;
    167147
    168   //OutputTestMolecule(TestMolecule1, "testmolecule1.xyz");
     148  //OutputTestMolecule1(TestMolecule1, "testmolecule1.xyz");
     149
     150  // create TestMolecule2 as copy
     151  TestMolecule2 = TestMolecule1->CopyMolecule();
     152  molecules->insert(TestMolecule2);
     153  CPPUNIT_ASSERT_EQUAL( TestMolecule2->getAtomCount(), 3 );
    169154
    170155  cout << "Case 1: offset of (3,0,0), hence angles are (104.5, 0, 75.5, 180) < 30." << endl;
    171156  Translator  = Vector(3,0,0);
    172   TestMolecule2->Translate(&Translator);
     157  TestMolecule1->Translate(&Translator);
    173158  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    174159  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, oxygen, NULL) );
    175   //OutputTestMolecule(TestMolecule2, "testmolecule2-1.xyz");
     160  //OutputTestMolecule1(TestMolecule1, "testmolecule2-1.xyz");
    176161  Translator = Vector(-3,0,0);
    177   TestMolecule2->Translate(&Translator);
     162  TestMolecule1->Translate(&Translator);
    178163
    179164  cout << "Case 2: offset of (0,3,0), hence angle are (14.5, 165.5, 90) < 30 (only three, because other 90 is missing due to first H01 only fulfilling H-bond criteria)." << endl;
    180165  Translator = Vector(0,3,0);
    181   TestMolecule2->Translate(&Translator);
     166  TestMolecule1->Translate(&Translator);
    182167  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    183   //OutputTestMolecule(TestMolecule2, "testmolecule2-2.xyz");
     168  //OutputTestMolecule1(TestMolecule1, "testmolecule2-2.xyz");
    184169  Translator = Vector(0,-3,0);
    185   TestMolecule2->Translate(&Translator);
     170  TestMolecule1->Translate(&Translator);
    186171
    187172  cout << "Case 3: offset of (0,-3,0) and mirror, hence angle are (165.5, 90, 165.5, 90) > 30." << endl;
    188173  Translator = Vector(0,-3,0);
    189   TestMolecule2->Scale((const double ** const)&mirror);
    190   TestMolecule2->Translate(&Translator);
    191   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    192   //OutputTestMolecule(TestMolecule2, "testmolecule2-3.xyz");
     174  TestMolecule1->Scale((const double ** const)&mirror);
     175  TestMolecule1->Translate(&Translator);
     176  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
     177  //OutputTestMolecule1(TestMolecule1, "testmolecule2-3.xyz");
    193178  Translator = Vector(0,3,0);
    194   TestMolecule2->Translate(&Translator);
    195   TestMolecule2->Scale((const double ** const)&mirror);
     179  TestMolecule1->Translate(&Translator);
     180  TestMolecule1->Scale((const double ** const)&mirror);
    196181
    197182  cout << "Case 4: offset of (2,1,0), hence angle are (78, 26.6, 102, 153.4) < 30." << endl;
    198183  Translator = Vector(2,1,0);
    199   TestMolecule2->Translate(&Translator);
     184  TestMolecule1->Translate(&Translator);
    200185  CPPUNIT_ASSERT_EQUAL( 1 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    201   //OutputTestMolecule(TestMolecule2, "testmolecule2-4.xyz");
     186  //OutputTestMolecule1(TestMolecule1, "testmolecule2-4.xyz");
    202187  Translator = Vector(-2,-1,0);
    203   TestMolecule2->Translate(&Translator);
     188  TestMolecule1->Translate(&Translator);
    204189
    205190  cout << "Case 5: offset of (0,0,3), hence angle are (90, 90, 90, 90) > 30." << endl;
    206191  Translator = Vector(0,0,3);
    207   TestMolecule2->Translate(&Translator);
    208   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    209   //OutputTestMolecule(TestMolecule2, "testmolecule2-5.xyz");
     192  TestMolecule1->Translate(&Translator);
     193  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
     194  //OutputTestMolecule1(TestMolecule1, "testmolecule2-5.xyz");
    210195  Translator = Vector(0,0,-3);
    211   TestMolecule2->Translate(&Translator);
     196  TestMolecule1->Translate(&Translator);
    212197
    213198  cout << "Case 6: offset of (-3,0,0) and mirror, hence angle are (75.5, 180, 104.5, 180) > 30." << endl;
    214199  Translator = Vector(-3,0,0);
    215   TestMolecule2->Scale((const double ** const)&mirror);
    216   TestMolecule2->Translate(&Translator);
    217   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    218   //OutputTestMolecule(TestMolecule2, "testmolecule2-6.xyz");
     200  TestMolecule1->Scale((const double ** const)&mirror);
     201  TestMolecule1->Translate(&Translator);
     202  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
     203  //OutputTestMolecule1(TestMolecule1, "testmolecule2-6.xyz");
    219204  Translator = Vector(3,0,0);
    220   TestMolecule2->Translate(&Translator);
    221   TestMolecule2->Scale((const double ** const)&mirror);
     205  TestMolecule1->Translate(&Translator);
     206  TestMolecule1->Scale((const double ** const)&mirror);
    222207
    223208  cout << "Case 7: offset of (3,0,0) and mirror, hence angles are (104.5, 0, 104.5, 0) < 30, but interfering hydrogens." << endl;
    224209  Translator = Vector(3,0,0);
    225   TestMolecule2->Scale((const double ** const)&mirror);
    226   TestMolecule2->Translate(&Translator);
    227   CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    228   //OutputTestMolecule(TestMolecule2, "testmolecule2-7.xyz");
     210  TestMolecule1->Scale((const double ** const)&mirror);
     211  TestMolecule1->Translate(&Translator);
     212  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
     213  //OutputTestMolecule1(TestMolecule1, "testmolecule2-7.xyz");
    229214  Translator = Vector(-3,0,0);
    230   TestMolecule2->Translate(&Translator);
    231   TestMolecule2->Scale((const double ** const)&mirror);
     215  TestMolecule1->Translate(&Translator);
     216  TestMolecule1->Scale((const double ** const)&mirror);
    232217
    233218  cout << "Case 8: offset of (0,3,0), hence angle are (14.5, 90, 14.5, 90) < 30, but interfering hydrogens." << endl;
    234219  Translator = Vector(0,3,0);
    235   TestMolecule2->Scale((const double ** const)&mirror);
    236   TestMolecule2->Translate(&Translator);
    237   //OutputTestMolecule(TestMolecule2, "testmolecule2-8.xyz");
     220  TestMolecule1->Scale((const double ** const)&mirror);
     221  TestMolecule1->Translate(&Translator);
     222  //OutputTestMolecule1(TestMolecule1, "testmolecule2-8.xyz");
    238223  CPPUNIT_ASSERT_EQUAL( 0 , CountHydrogenBridgeBonds(molecules, NULL, NULL) );
    239224  Translator = Vector(0,-3,0);
    240   TestMolecule2->Translate(&Translator);
    241   TestMolecule2->Scale((const double ** const)&mirror);
     225  TestMolecule1->Translate(&Translator);
     226  TestMolecule1->Scale((const double ** const)&mirror);
    242227
    243228  delete[](mirror);

  • src/Graph/BondGraph.cpp

    rcbcda6 r826e8c  
    9898    int secondZ) const
    9999{
     100  double return_length;
     101  if ((firstZ < 0) || (firstZ >= (int)BondLengthMatrix->Matrix[0].size()))
     102    return -1.;
     103  if ((secondZ < 0) || (secondZ >= (int)BondLengthMatrix->Matrix[0][firstZ].size()))
     104    return -1.;
    100105  if (BondLengthMatrix == NULL) {
    101     std::cout << "-1." << std::endl;
    102     return( -1. );
     106    return_length = -1.;
    103107  } else {
    104     std::cout << BondLengthMatrix->Matrix[0][firstZ][secondZ] << std::endl;
    105     return (BondLengthMatrix->Matrix[0][firstZ][secondZ]);
    106   }
     108    return_length = BondLengthMatrix->Matrix[0][firstZ][secondZ];
     109  }
     110  LOG(4, "INFO: Request for length between " << firstZ << " and "
     111      << secondZ << ": " << return_length << ".");
     112  return return_length;
    107113}
    108114
    … …  
    168174}
    169175
    170 void BondGraph::CreateAdjacency(LinkedCell_deprecated &LC) const
    171 {
    172   atom *Walker = NULL;
    173   atom *OtherWalker = NULL;
    174   int n[NDIM];
    175   Box &domain = World::getInstance().getDomain();
    176   size_t CurrentTime = WorldTime::getTime();
    177 
    178   unsigned int BondCount = 0;
    179   // 3a. go through every cell
    180   LOG(3, "INFO: Celling ... ");
    181   for (LC.n[0] = 0; LC.n[0] < LC.N[0]; LC.n[0]++)
    182     for (LC.n[1] = 0; LC.n[1] < LC.N[1]; LC.n[1]++)
    183       for (LC.n[2] = 0; LC.n[2] < LC.N[2]; LC.n[2]++) {
    184         const TesselPointSTLList *List = LC.GetCurrentCell();
    185         LOG(2, "Current cell is " << LC.n[0] << ", " << LC.n[1] << ", " << LC.n[2] << " with No. " << LC.index << " containing " << List->size() << " points.");
    186         if (List != NULL) {
    187           for (TesselPointSTLList::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    188             Walker = dynamic_cast<atom*>(*Runner);
    189             ASSERT(Walker != NULL,
    190                 "BondGraph::CreateAdjacency() - Tesselpoint that was not an atom retrieved from LinkedNode");
    191             LOG(2, "INFO: Current Atom is " << *Walker << ".");
    192             // 3c. check for possible bond between each atom in this and every one in the 27 cells
    193             for (n[0] = -1; n[0] <= 1; n[0]++)
    194               for (n[1] = -1; n[1] <= 1; n[1]++)
    195                 for (n[2] = -1; n[2] <= 1; n[2]++) {
    196                   const TesselPointSTLList *OtherList = LC.GetRelativeToCurrentCell(n);
    197                   if (OtherList != NULL) {
    198                     LOG(3, "INFO: Current relative cell is " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << LC.index << " containing " << OtherList->size() << " points.");
    199                     for (TesselPointSTLList::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
    200                       if ((*OtherRunner) > Walker) {  // just to not add bonds from both sides
    201                         OtherWalker = dynamic_cast<atom*>(*OtherRunner);
    202                         ASSERT(OtherWalker != NULL,
    203                             "BondGraph::CreateAdjacency() - TesselPoint that was not an atom retrieved from LinkedNode");
    204                         LOG(3, "INFO: Current other atom is " << *OtherWalker << ".");
    205                         if (OtherWalker->father > Walker->father ) { // just to not add bonds from both sides
    206                           const range<double> MinMaxDistanceSquared(
    207                               getMinMaxDistanceSquared(Walker, OtherWalker));
    208                           const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
    209                           LOG(3, "INFO: Checking squared distance " << distance << " against typical bond length of " << MinMaxDistanceSquared << ".");
    210                           const bool status = MinMaxDistanceSquared.isInRange(distance);
    211                           if (status) { // create bond if distance is smaller
    212                             LOG(1, "ACCEPT: Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << ".");
    213                             //const bond * Binder =
    214                               Walker->father->addBond(CurrentTime, OtherWalker->father);
    215                             BondCount++;
    216                           } else {
    217                             LOG(2, "REJECT: Squared distance "
    218                                 << distance << " is out of squared covalent bounds "
    219                                 << MinMaxDistanceSquared << ".");
    220                           }
    221                         } else {
    222                           LOG(5, "REJECT: Not Adding: Wrong order of father's.");
    223                         }
    224                       } else {
    225                         LOG(5, "REJECT: Not Adding: Wrong order.");
    226                       }
    227                     }
    228                   }
    229                 }
    230           }
    231         }
    232       }
    233     LOG(1, "I detected " << BondCount << " bonds in the molecule.");
     176LinkedCell::LinkedCell_View BondGraph::getLinkedCell(const double max_distance) const
     177{
     178  return World::getInstance().getLinkedCell(max_distance);
     179}
     180
     181Box &BondGraph::getDomain() const
     182{
     183  return World::getInstance().getDomain();
     184}
     185
     186unsigned int BondGraph::getTime() const
     187{
     188  return WorldTime::getTime();
    234189}
    235190

  • src/Graph/BondGraph.hpp

    rcbcda6 r826e8c  
    2424#include "Atom/AtomSet.hpp"
    2525#include "Bond/bond.hpp"
     26#include "Box.hpp"
    2627#include "CodePatterns/Assert.hpp"
    2728#include "CodePatterns/Log.hpp"
    … …  
    2930#include "Element/element.hpp"
    3031#include "Fragmentation/MatrixContainer.hpp"
    31 #include "LinkedCell/linkedcell.hpp"
     32#include "Helpers/defs.hpp"
     33#include "LinkedCell/LinkedCell_View.hpp"
    3234#include "LinkedCell/IPointCloud.hpp"
    3335#include "LinkedCell/PointCloudAdaptor.hpp"
    … …  
    107109  }
    108110
     111  /** Determines the maximum of all element::CovalentRadius for elements present in \a &Set.
     112   *
     113   * I.e. the function returns a sensible cutoff criteria for bond recognition,
     114   * e.g. to be used for LinkedCell_deprecated or others.
     115   *
     116   * \param &Set AtomSetMixin with all particles to consider
     117   */
     118  template <class container_type,
     119            class iterator_type,
     120            class const_iterator_type>
     121  double getMaxPossibleBondDistance(
     122      const element * const Walker,
     123      const AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set) const
     124  {
     125    double max_distance = 0.;
     126    // get all elements
     127    std::set< const element *> PresentElements;
     128    for(const_iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
     129      PresentElements.insert( (*AtomRunner)->getType() );
     130    }
     131    // create all element combinations
     132    for (std::set< const element *>::const_iterator iter = PresentElements.begin();
     133        iter != PresentElements.end();
     134        ++iter) {
     135      const range<double> MinMaxDistance(getMinMaxDistance((*iter),Walker));
     136      if (MinMaxDistance.last > max_distance)
     137        max_distance = MinMaxDistance.last;
     138    }
     139    return max_distance;
     140  }
     141
    109142  /** Returns bond criterion for given pair based on a bond length matrix.
    110143   * This calls element-version of getMinMaxDistance().
    … …  
    128161      const BondedParticle * const OtherWalker) const;
    129162
    130   /** Creates the adjacency list for a given \a Range of iterable atoms.
    131    *
    132    * @param Set Range with begin and end iterator
     163  /** Creates an adjacency list of the molecule.
     164   * Generally, we use the CSD approach to bond recognition, that is the the distance
     165   * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
     166   * a threshold t = 0.4 Angstroem.
     167   * To make it O(N log N) the function uses the linked-cell technique as follows:
     168   * The procedure is step-wise:
     169   *  -# Remove every bond in list
     170   *  -# go through every atom in given \a set, check the atoms therein against all possible bond partners, add bond if true
     171   *  -# correct the bond degree iteratively (single->double->triple bond)
     172   *  -# finally print the bond list to \a *out if desired
     173   * \param &set Container with all atoms to create adjacency for
    133174   */
    134175  template <class container_type,
    … …  
    145186    if (counter > 1) {
    146187      LOG(1, "STATUS: Setting max bond distance.");
    147       const double max_distance = getMaxPossibleBondDistance(Set);
     188      LinkedCell::LinkedCell_View LC = getLinkedCell(getMaxPossibleBondDistance(Set));
    148189
    149190      LOG(1, "STATUS: Creating LinkedCell structure for given " << counter << " atoms.");
    150       PointCloudAdaptor< AtomSetMixin<container_type,iterator_type> > cloud(&Set, "SetOfAtoms");
    151       LinkedCell_deprecated *LC = new LinkedCell_deprecated(cloud, max_distance);
    152 
    153       CreateAdjacency(*LC);
    154       delete (LC);
     191
     192      Box &domain = getDomain();
     193      unsigned int CurrentTime = getTime();
     194
     195      unsigned int BondCount = 0;
     196      // go through every atom in the set (observed cause we change its bonds)
     197      for(typename AtomSetMixin<container_type,iterator_type,const_iterator_type>::iterator iter = Set.begin();
     198          iter != Set.end(); ++iter) {
     199        const atom * const Walker = dynamic_cast<const atom *>(*iter);
     200        ASSERT(Walker != NULL,
     201            "BondGraph::CreateAdjacency() - TesselPoint "
     202            +(*iter)->getName()+" that was not an atom retrieved from given set");
     203        LOG(2, "INFO: Current Atom is " << *Walker << ".");
     204
     205        // obtain all possible neighbors
     206        LinkedCell::LinkedList ListOfNeighbors = LC.getAllNeighbors(
     207            getMaxPossibleBondDistance(Walker->getType(), Set),
     208            Walker->getPosition());
     209        if (!ListOfNeighbors.empty()) {
     210          // we have some possible candidates, go through each
     211          for (LinkedCell::LinkedList::const_iterator neighboriter = ListOfNeighbors.begin();
     212              neighboriter != ListOfNeighbors.end();
     213              ++neighboriter) {
     214            if ((*neighboriter) > Walker) {  // just to not add bonds from both sides
     215              const atom * const OtherWalker = dynamic_cast<const atom *>(*neighboriter);
     216              ASSERT(OtherWalker != NULL,
     217                  "BondGraph::CreateAdjacency() - TesselPoint "
     218                  +(*neighboriter)->getName()+" that was not an atom retrieved from LinkedList");
     219              LOG(3, "INFO: Current other atom is " << *OtherWalker << ".");
     220
     221              const range<double> MinMaxDistanceSquared(
     222                  getMinMaxDistanceSquared(Walker, OtherWalker));
     223              const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
     224              LOG(3, "INFO: Checking squared distance " << distance << " against typical bond length of " << MinMaxDistanceSquared << ".");
     225              const bool status = MinMaxDistanceSquared.isInRange(distance);
     226              if (status) { // create bond if distance is smaller
     227                LOG(1, "ACCEPT: Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << ".");
     228                // directly use iter to avoid const_cast'ing Walker, too
     229                //const bond * Binder =
     230                const_cast<atom *>(Walker)->addBond(CurrentTime, const_cast<atom *>(OtherWalker));
     231                ++BondCount;
     232              } else {
     233                LOG(2, "REJECT: Squared distance "
     234                    << distance << " is out of squared covalent bounds "
     235                    << MinMaxDistanceSquared << ".");
     236              }
     237
     238            } else {
     239              LOG(5, "REJECT: Not Adding: Wrong order.");
     240            }
     241          }
     242        }
     243      }
     244      LOG(1, "I detected " << BondCount << " bonds in the molecule.");
    155245
    156246      // correct bond degree by comparing valence and bond degree
    … …  
    242332    LOG(1, "STATUS: "<< bondcounter << " bonds have been parsed.");
    243333  }
    244 
    245   /** Creates an adjacency list of the molecule.
    246    * Generally, we use the CSD approach to bond recognition, that is the the distance
    247    * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
    248    * a threshold t = 0.4 Angstroem.
    249    * To make it O(N log N) the function uses the linked-cell technique as follows:
    250    * The procedure is step-wise:
    251    *  -# Remove every bond in list
    252    *  -# Count the atoms in the molecule with CountAtoms()
    253    *  -# partition cell into smaller linked cells of size \a bonddistance
    254    *  -# put each atom into its corresponding cell
    255    *  -# go through every cell, check the atoms therein against all possible bond partners in the 27 adjacent cells, add bond if true
    256    *  -# correct the bond degree iteratively (single->double->triple bond)
    257    *  -# finally print the bond list to \a *out if desired
    258    * \param &LC Linked Cell Container with all atoms
    259    */
    260   void CreateAdjacency(LinkedCell_deprecated &LC) const;
    261334
    262335  /** Removes all bonds within the given set of iterable atoms.
    … …  
    325398
    326399private:
     400  /** Convenience function to place access to World::getLinkedCell() into source module.
     401   *
     402   * @return ref to LinkedCell_View
     403   */
     404  LinkedCell::LinkedCell_View getLinkedCell(const double max_distance) const;
     405
     406  /** Convenience function to place access to World::getDomain() into source module.
     407   *
     408   * @return ref to Box
     409   */
     410  Box &getDomain() const;
     411
     412  /** Convenience function to place access to WorldTime::getTime() into source module.
     413   *
     414   * @return current time step
     415   */
     416  unsigned int getTime() const;
    327417
    328418  /** Returns the BondLengthMatrix entry for a given index pair.
Note: See TracChangeset for help on using the changeset viewer.