Changeset ea7176 for src/molecule.cpp

Timestamp:

Mar 19, 2010, 1:29:01 PM (15 years ago)

Author:

Tillmann Crueger <crueger@…>

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:

e0b6fd

Parents:

80c63d

Message:

FIX: Made AtomCount a Cacheable so that the number of atoms in a molecule will always be correct

All unittests working.
All Complete testcases fail.

File:

• src/molecule.cpp (modified) (15 diffs)

src/molecule.cpp

@@ -32 +32 @@
  */
 molecule::molecule(const periodentafel * const teil) : elemente(teil),
-  first(new bond(0, 0, 1, -1)), last(new bond(0, 0, 1, -1)), MDSteps(0), AtomCount(0),
+  first(new bond(0, 0, 1, -1)), last(new bond(0, 0, 1, -1)), MDSteps(0),
   BondCount(0), ElementCount(0), NoNonHydrogen(0), NoNonBonds(0), NoCyclicBonds(0), BondDistance(0.), 
   ActiveFlag(false), IndexNr(-1),
   formula(this,boost::bind(&molecule::calcFormula,this)),
-  last_atom(0),
-  InternalPointer(begin())
+  AtomCount(this,boost::bind(&molecule::doCountAtoms,this)), last_atom(0),  InternalPointer(begin())
 {
   // init bond chain list
@@ -72 +71 @@
 const std::string molecule::getName(){
   return std::string(name);
+}
+
+int molecule::getAtomCount() const{
+  return *AtomCount;
 }
 
@@ -180 +183 @@
   if (pointer != NULL) {
     pointer->sort = &pointer->nr;
-    AtomCount++;
     if (pointer->type != NULL) {
       if (ElementsInMolecule[pointer->type->Z] == 0)
@@ -214 +216 @@
     if ((pointer->type != NULL) && (pointer->type->Z != 1))
       NoNonHydrogen++;
-    AtomCount++;
     retval=walker;
   }
@@ -611 +612 @@
 
   // copy values
-  copy->CountAtoms();
   copy->CountElements();
   if (first->next != last) {  // if adjaceny list is present
@@ -728 +728 @@
 bool molecule::RemoveAtom(atom *pointer)
 {
+  OBSERVE;
   if (ElementsInMolecule[pointer->type->Z] != 0)  { // this would indicate an error
     ElementsInMolecule[pointer->type->Z]--;  // decrease number of atom of this element
-    AtomCount--;
   } else
     eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
@@ -909 +909 @@
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
     for (int step=0;step<MDSteps;step++) {
-      *output << AtomCount << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
+      *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
       ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step );
     }
@@ -926 +926 @@
   if (output != NULL) {
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
-    *output << AtomCount << "\n\tCreated by molecuilder on " << ctime(&now);
+    *output << getAtomCount() << "\n\tCreated by molecuilder on " << ctime(&now);
     ActOnAllAtoms( &atom::OutputXYZLine, output );
     return true;
@@ -936 +936 @@
  * \param *out output stream for debugging
  */
-void molecule::CountAtoms()
-{
-  int i = size();
-  if ((AtomCount == 0) || (i != AtomCount)) {
-    cout << "!!!!!!!!! Counting needed" << endl;
-    Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
-    AtomCount = i;
-
-    // count NonHydrogen atoms and give each atom a unique name
-    if (AtomCount != 0) {
-      i=0;
-      NoNonHydrogen = 0;
-      for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-        (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
-        if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
-          NoNonHydrogen++;
-        Free(&(*iter)->Name);
-        (*iter)->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
-        sprintf((*iter)->Name, "%2s%02d", (*iter)->type->symbol, (*iter)->nr+1);
-        Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->Name << "." << endl;
-        i++;
-      }
-    } else
-      Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
-  }
+int molecule::doCountAtoms()
+{
+  int res = size();
+  int i = 0;
+  NoNonHydrogen = 0;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
+    if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
+      NoNonHydrogen++;
+    Free(&(*iter)->Name);
+    (*iter)->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
+    sprintf((*iter)->Name, "%2s%02d", (*iter)->type->symbol, (*iter)->nr+1);
+    Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->Name << "." << endl;
+    i++;
+  }
+  return res;
 };
 
@@ -1022 +1013 @@
   /// first count both their atoms and elements and update lists thereby ...
   //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
-  CountAtoms();
-  OtherMolecule->CountAtoms();
   CountElements();
   OtherMolecule->CountElements();
@@ -1030 +1019 @@
   /// -# AtomCount
   if (result) {
-    if (AtomCount != OtherMolecule->AtomCount) {
-      Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
+    if (getAtomCount() != OtherMolecule->getAtomCount()) {
+      Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
       result = false;
-    } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
+    } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
   }
   /// -# ElementCount
@@ -1073 +1062 @@
   if (result) {
     Log() << Verbose(5) << "Calculating distances" << endl;
-    Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
-    OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
+    Distances = Calloc<double>(getAtomCount(), "molecule::IsEqualToWithinThreshold: Distances");
+    OtherDistances = Calloc<double>(getAtomCount(), "molecule::IsEqualToWithinThreshold: OtherDistances");
     SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
     SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
@@ -1080 +1069 @@
     /// ... sort each list (using heapsort (o(N log N)) from GSL)
     Log() << Verbose(5) << "Sorting distances" << endl;
-    PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
-    OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
-    gsl_heapsort_index (PermMap, Distances, AtomCount, sizeof(double), CompareDoubles);
-    gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
-    PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
+    PermMap = Calloc<size_t>(getAtomCount(), "molecule::IsEqualToWithinThreshold: *PermMap");
+    OtherPermMap = Calloc<size_t>(getAtomCount(), "molecule::IsEqualToWithinThreshold: *OtherPermMap");
+    gsl_heapsort_index (PermMap, Distances, getAtomCount(), sizeof(double), CompareDoubles);
+    gsl_heapsort_index (OtherPermMap, OtherDistances, getAtomCount(), sizeof(double), CompareDoubles);
+    PermutationMap = Calloc<int>(getAtomCount(), "molecule::IsEqualToWithinThreshold: *PermutationMap");
     Log() << Verbose(5) << "Combining Permutation Maps" << endl;
-    for(int i=AtomCount;i--;)
+    for(int i=getAtomCount();i--;)
       PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
 
@@ -1092 +1081 @@
     Log() << Verbose(4) << "Comparing distances" << endl;
     flag = 0;
-    for (int i=0;i<AtomCount;i++) {
+    for (int i=0;i<getAtomCount();i++) {
       Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
       if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
@@ -1129 +1118 @@
 {
   Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
-  int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
-  for (int i=AtomCount;i--;)
+  int *AtomicMap = Malloc<int>(getAtomCount(), "molecule::GetAtomicMap: *AtomicMap");
+  for (int i=getAtomCount();i--;)
     AtomicMap[i] = -1;
   if (OtherMolecule == this) {  // same molecule
-    for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
+    for (int i=getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence
       AtomicMap[i] = i;
     Log() << Verbose(4) << "Map is trivial." << endl;