Changeset 274d45 for src

Timestamp:

Jun 19, 2010, 3:54:47 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:

bf0edf

Parents:

ad2b411

git-author:

Frederik Heber <heber@…> (06/19/10 12:35:01)

git-committer:

Frederik Heber <heber@…> (06/19/10 15:54:47)

Message:

FIX: Atoms were stored not in the sequence they were loaded.

1. The main problem is molecule::atomSet which is a set<atom *>, i.e. atoms are sorted by their appearance in memory. As memory need not be allocated sequentially, this gives rise to extreme arbitririty which is not desired. Instead the atoms should be stored in the sequence they were loaded/created. The solution is as follows:

• config::SaveAll()
  • uses SimpleMultiMerge() instead of SimpleMultiAdd(). I.e. no atoms are copied (with new ids) but are relinked into a single molecule.
  • also, in case of only a single molecule in the list, we don't merge anything to avoid mol id changes.
  • MaxMinSteps were set according to MaxPsiDouble overriding correct results from molecule::CalculateOrbitals()
• molecule::atomSet is now a list<atom *>
• molecule::atomIds is a new set<atomId_t> (atomIdSet) which controls that (global) ids remain unique in the no more Atomset's set (but list)
• molecule::erase() erases also in atomIds
• molecule::insert() checks whether id is present by atomIds
• molecule::find() as std::list does not have a find, we just go through the list until the object is found (or not), this may be speeded up by another internal list.
• molecule::InternalPointer made lots of confusion as virtual function GoToFirst() is const, hence begin() (needed therein) returns const_iterator, which then cannot be simply re-cast into an iterator: We make it a pointer, reinterpret_cast the pointer and reference it back. Although InternalPointer is mutable, the compiler cannot use the non-const function begin() (it cannot be made const, as overloading is not allowed). (this is noted in the code also extensively.)
• molecule::containsAtom() does not use count but the new find, as it returns only boolean anyway.
• rewrote MoleculeListClass::SimpleMerge() to get rid of the extra iterator, as we remove all atoms in the end anyway.
• FIX: MoleculeListClass::SimpleMultiMerge() - the created mol was not inserted into the moleculelist in the end, although it is the only one to remain.

2. All other databases had missing headers with respect to those stored in elements_db.cpp. Hence, valence of hydrogen was not parsed and this caused several failures in CalculateOrbitals() (Psi numbers and MaxMinSetp in pcp conf file).

3. Subsequenytly, various test cases (12, 21, 30, 31, 36-38, 39) were broken. This had two reasons:

• Seemingly, CalculateOrbitals() was broken before hence MaxMinStep (pcp config) and MaxPsiDouble/PsiMaxNn[Up|Down] were always 0. (10-21,30-31,39)
• As the order is now correct, fixes from commits c9217161ec2a5d5db508557fe98a32068461f45b and 22a6da8380911571debebd69444d2615450bbbd8 were obselete and have been reverted (order of the Ion?_Type...): Molecules/6 (30), Molecules/7 (31), Filling/1 (39)
• Due to different ordering, Tesselation/3 (38) had completely different .dat file (though same tesselation)
• r3d had small differences, mostly order or epsilon (0 not being 0 but ..-e16), hence diffing was deactivated, as r3d is deprecated anyway (since vmd can render triangles as well and better).

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

Location:

src

Files:

• Hbondangle.db (modified) (1 diff)
• Hbonddistance.db (modified) (1 diff)
• config.cpp (modified) (3 diffs)
• molecule.cpp (modified) (7 diffs)
• molecule.hpp (modified) (2 diffs)
• molecule_pointcloud.cpp (modified) (1 diff)
• moleculelist.cpp (modified) (2 diffs)
• orbitals.db (modified) (1 diff)
• periodentafel.cpp (modified) (2 diffs)
• valence.db (modified) (1 diff)

src/Hbondangle.db

@@ +1 @@
+# atomicnumber angles for single, double and triple bond (-1 no angle)
 1       180     -1      -1
 5       180     131.0   109.2

src/Hbonddistance.db

@@ +1 @@
+#atomic number bond distances for single, double and triple bond (-1 no bond)
 1       0.74    -1      -1
 2       0.77429209      -1      -1

src/config.cpp

@@ -1784 +1784 @@
   char filename[MAXSTRINGSIZE];
   ofstream output;
-  molecule *mol = World::getInstance().createMolecule();
-  mol->SetNameFromFilename(ConfigFileName);
+  molecule *mol = NULL;
 
   if (!strcmp(configpath, GetDefaultPath())) {
@@ -1816 +1815 @@
   // translate each to its center and merge all molecules in MoleculeListClass into this molecule
   int N = molecules->ListOfMolecules.size();
-  int *src = new int[N];
-  N=0;
-  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
-    src[N++] = (*ListRunner)->IndexNr;
-    (*ListRunner)->Translate(&(*ListRunner)->Center);
-  }
-  molecules->SimpleMultiAdd(mol, src, N);
-  delete[](src);
-
-  // ... and translate back
-  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
-    (*ListRunner)->Center.Scale(-1.);
-    (*ListRunner)->Translate(&(*ListRunner)->Center);
-    (*ListRunner)->Center.Scale(-1.);
+  if (N != 1) { // don't do anything in case of only one molecule (shifts mol ids otherwise)
+    int *src = new int[N];
+    N=0;
+    for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      src[N++] = (*ListRunner)->IndexNr;
+      (*ListRunner)->Translate(&(*ListRunner)->Center);
+    }
+    mol = World::getInstance().createMolecule();
+    mol->SetNameFromFilename(ConfigFileName);
+    molecules->SimpleMultiMerge(mol, src, N);
+    mol->doCountAtoms();
+    mol->CountElements();
+    mol->CalculateOrbitals(*this);
+    delete[](src);
+  } else {
+    if (!molecules->ListOfMolecules.empty()) {
+      mol = *(molecules->ListOfMolecules.begin());
+      mol->doCountAtoms();
+      mol->CalculateOrbitals(*this);
+    } else {
+      DoeLog(1) && (eLog() << Verbose(1) << "There are no molecules to save!" << endl);
+    }
   }
 
   Log() << Verbose(0) << "Storing configuration ... " << endl;
   // get correct valence orbitals
-  mol->CalculateOrbitals(*this);
-  InitMaxMinStopStep = MaxMinStopStep = MaxPsiDouble;
   if (ConfigFileName != NULL) { // test the file name
     strcpy(filename, ConfigFileName);
@@ -1895 +1900 @@
   }
 
-  World::getInstance().destroyMolecule(mol);
+  // don't destroy molecule as it contains all our atoms
+  //World::getInstance().destroyMolecule(mol);
 };
 

src/molecule.cpp

@@ -42 +42 @@
   NoCyclicBonds(0), BondDistance(0.),  ActiveFlag(false), IndexNr(-1),
   formula(this,boost::bind(&molecule::calcFormula,this),"formula"),
-  AtomCount(this,boost::bind(&molecule::doCountAtoms,this),"AtomCount"), last_atom(0),  InternalPointer(begin())
+  AtomCount(this,boost::bind(&molecule::doCountAtoms,this),"AtomCount"), last_atom(0),  InternalPointer(atoms.begin())
 {
 
@@ -146 +146 @@
   iter--;
   atom* atom = *loc;
-  atoms.erase( loc );
+  atomIds.erase( atom->getId() );
+  atoms.remove( atom );
   atom->removeFromMolecule();
   return iter;
@@ -156 +157 @@
   molecule::const_iterator iter = find(key);
   if (iter != end()){
-    atoms.erase( iter++ );
+    atomIds.erase( key->getId() );
+    atoms.remove( key );
     key->removeFromMolecule();
   }
@@ -164 +166 @@
 molecule::const_iterator molecule::find ( atom * key ) const
 {
-  return atoms.find( key );
+  molecule::const_iterator iter;
+  for (molecule::const_iterator Runner = begin(); Runner != end(); ++Runner) {
+    if (*Runner == key)
+      return molecule::const_iterator(Runner);
+  }
+  return molecule::const_iterator(atoms.end());
 }
 
 pair<molecule::iterator,bool> molecule::insert ( atom * const key )
 {
-  pair<atomSet::iterator,bool> res = atoms.insert(key);
-  return pair<iterator,bool>(iterator(res.first,this),res.second);
+  pair<atomIdSet::iterator,bool> res = atomIds.insert(key->getId());
+  if (res.second) { // push atom if went well
+    atoms.push_back(key);
+    return pair<iterator,bool>(molecule::iterator(--end()),res.second);
+  } else {
+    return pair<iterator,bool>(molecule::iterator(end()),res.second);
+  }
 }
 
 bool molecule::containsAtom(atom* key){
-  return atoms.count(key);
+  return (find(key) != end());
 }
 
@@ -790 +802 @@
   for (molecule::iterator iter = begin(); !empty(); iter = begin())
       erase(iter);
+  return empty();
 };
 
@@ -1010 +1023 @@
   configuration.MaxPsiDouble /= 2;
   configuration.PsiType = (configuration.PsiMaxNoDown == configuration.PsiMaxNoUp) ? 0 : 1;
-  if ((configuration.PsiType == 1) && (configuration.ProcPEPsi < 2)) {
+  if ((configuration.PsiType == 1) && (configuration.ProcPEPsi < 2) && ((configuration.PsiMaxNoDown != 1) || (configuration.PsiMaxNoUp != 0))) {
     configuration.ProcPEGamma /= 2;
     configuration.ProcPEPsi *= 2;
@@ -1017 +1030 @@
     configuration.ProcPEPsi = 1;
   }
-  configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.MaxPsiDouble;
+  cout << configuration.PsiMaxNoDown << ">" << configuration.PsiMaxNoUp << endl;
+  if (configuration.PsiMaxNoDown > configuration.PsiMaxNoUp) {
+    configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.PsiMaxNoDown;
+    cout << configuration.PsiMaxNoDown << " " << configuration.InitMaxMinStopStep << endl;
+  } else {
+    configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.PsiMaxNoUp;
+    cout << configuration.PsiMaxNoUp << " " << configuration.InitMaxMinStopStep << endl;
+  }
 };
 

src/molecule.hpp

@@ -93 +93 @@
 
   public:
-    typedef std::set<atom*> atomSet;
+    typedef std::list<atom*> atomSet;
+    typedef std::set<atomId_t> atomIdSet;
     typedef ObservedIterator<atomSet> iterator;
     typedef atomSet::const_iterator const_iterator;
@@ -121 +122 @@
     Cacheable<int>    AtomCount;
     moleculeId_t id;
-    atomSet atoms; //<!set of atoms
+    atomSet atoms; //<!list of atoms
+    atomIdSet atomIds; //<!set of atomic ids to check uniqueness of atoms
   protected:
     //void CountAtoms();

src/molecule_pointcloud.cpp

@@ -57 +57 @@
 void molecule::GoToFirst() const
 {
-  InternalPointer = atoms.begin();
+  // evil hack necessary because
+  // -# although InternalPointer is mutable
+  // -# only const_iterator begin() is called due to const in the function declaration above
+  // -# and there is no cast from const_iterator to const iterator
+  atomSet::const_iterator test = begin();
+  InternalPointer = *(reinterpret_cast<atomSet::iterator *>(&test));
 };
 

src/moleculelist.cpp

@@ -212 +212 @@
 
   // put all molecules of src into mol
-  molecule::iterator runner;
-  for (molecule::iterator iter = srcmol->begin(); iter != srcmol->end(); ++iter) {
-    runner = iter++;
-    srcmol->UnlinkAtom((*runner));
-    mol->AddAtom((*runner));
+  for (molecule::iterator iter = srcmol->begin(); !srcmol->empty(); iter=srcmol->begin()) {
+    atom * const Walker = *iter;
+    srcmol->UnlinkAtom(Walker);
+    mol->AddAtom(Walker);
   }
 
@@ -259 +258 @@
     status = status && SimpleMerge(mol, srcmol);
   }
+  insert(mol);
   return status;
 };

src/orbitals.db

@@ +1 @@
+# atomicnumber numberoforbitals
 1       1
 2       0

src/periodentafel.cpp

@@ -373 +373 @@
       (*input) >> elements[Z]->Valence;
       (*input) >> ws;
-      //Log() << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->Valence << " valence electrons." << endl;
+      //Log() << Verbose(3) << "Element " << Z << " has " << FindElement(Z)->Valence << " valence electrons." << endl;
     }
     return true;
@@ -396 +396 @@
       (*input) >> elements[Z]->NoValenceOrbitals;
       (*input) >> ws;
-      //Log() << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->NoValenceOrbitals << " number of singly occupied valence orbitals." << endl;
+      //Log() << Verbose(3) << "Element " << Z << " has " << FindElement(Z)->NoValenceOrbitals << " number of singly occupied valence orbitals." << endl;
     }
     return true;

src/valence.db

@@ +1 @@
+#atomicnumber numberofvalenceorbitals
 1       0.10000000000000E+01
 2       0.20000000000000E+01