Changeset 66fd49 for src/boundary.cpp

Timestamp:

Feb 3, 2011, 9:51:18 AM (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, Candidate_v1.7.0, 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:

d6f886

Parents:

0b15bb

git-author:

Frederik Heber <heber@…> (12/30/10 20:52:17)

git-committer:

Frederik Heber <heber@…> (02/03/11 09:51:18)

Message:

Rewrite of FillVoidWithMoleculeAction.

FillVoidWithMoleculeAction:

• new parameter MinDistance and default value of 0.
• BUGFIX: filler is already created when parsing file, removed useless creation of it initially (also caused lots of confusion due to an "extra" molecule).
• Undo implemented, regression test inserted.
• Redo is somewhat hard to implement, as one would use performCall() if it only it would not retrieve its values from ValueStorage ...

FillVoidWithMolecule():

• filler is now the zeroth not the last molecule, marked by firstInsertion and firstInserter. Filler is removed if no molecules are filled.
• outsourced stuff into smaller functions
  • RandomizeMoleculePositions()
  • RemoveAtomsOutsideDomain()
  • isSpaceAroundPointVoid()
• removed FillIt to through every atom despite only CurrentPosition, indepedent of atom position, is checked.

TESTFIXES:

• Analysis/3: test.xyz changed because boundary is now 1.5 instead of 2.1 as 2.1 is not enough of molecules get filled in (and the filler already is).
• Analysis/3: tensid.data was actually lacking water at (0,0,0) which is after the rewrite present.

File:

• src/boundary.cpp (modified) (5 diffs)

src/boundary.cpp

@@ -927 +927 @@
 };
 
-/** Fills the empty space of the simulation box with water/
- * \param *out output stream for debugging
- * \param *TesselStruct contains tesselated surface
+/** Rotates given molecule \a Filling and moves its atoms according to given
+ *  \a RandomAtomDisplacement.
+ *
+ *  Note that for rotation to be sensible, the molecule should be centered at
+ *  the origin. This is not done here!
+ *
+ *  \param &Filling molecule whose atoms to displace
+ *  \param RandomAtomDisplacement magnitude of random displacement
+ *  \param &Rotations 3D rotation matrix (or unity if no rotation desired)
+ */
+void RandomizeMoleculePositions(
+    molecule *&Filling,
+    double RandomAtomDisplacement,
+    RealSpaceMatrix &Rotations
+    )
+{
+  Vector AtomTranslations;
+  for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ++miter) {
+    Vector temp = (*miter)->getPosition();
+    temp *= Rotations;
+    (*miter)->setPosition(temp);
+    // create atomic random translation vector ...
+    for (int i=0;i<NDIM;i++)
+      AtomTranslations[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+    (*miter)->setPosition((*miter)->getPosition() + AtomTranslations);
+  }
+}
+
+/** Removes all atoms of a molecule outside.
+ *
+ * If the molecule is empty, it is removed as well.
+ *
+ * @param Filling molecule whose atoms to check, removed if eventually left
+ *        empty.
+ */
+void RemoveAtomsOutsideDomain(molecule *&Filling)
+{
+  Box &Domain = World::getInstance().getDomain();
+  // check if all is still inside domain
+  for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ) {
+    // check whether each atom is inside box
+    if (!Domain.isInside((*miter)->getPosition())) {
+      atom *Walker = *miter;
+      ++miter;
+      World::getInstance().destroyAtom(Walker);
+    } else {
+      ++miter;
+    }
+  }
+  if (Filling->empty()) {
+    DoLog(0) && (Log() << Verbose(0) << "Removing molecule " << Filling->getName() << ", all atoms have been removed." << std::endl);
+    World::getInstance().destroyMolecule(Filling);
+  }
+}
+
+/** Checks whether there are no atoms inside a sphere around \a CurrentPosition
+ *  except those atoms present in \a *filler.
+ *  If filler is NULL, then we just call LinkedCell::GetPointsInsideSphere() and
+ *  check whether the return list is empty.
+ * @param *filler
+ * @param boundary
+ * @param CurrentPosition
+ */
+bool isSpaceAroundPointVoid(
+    LinkedCell *LC,
+    molecule *filler,
+    const double boundary,
+    Vector &CurrentPosition)
+{
+  size_t compareTo = 0;
+  LinkedCell::LinkedNodes* liste = LC->GetPointsInsideSphere(boundary == 0. ? MYEPSILON : boundary, &CurrentPosition);
+  if (filler != NULL) {
+    for (LinkedCell::LinkedNodes::const_iterator iter = liste->begin();
+        iter != liste->end();
+        ++iter) {
+      for (molecule::iterator miter = filler->begin();
+          miter != filler->end();
+          ++miter) {
+        if (*iter == *miter)
+          ++compareTo;
+      }
+    }
+  }
+  const bool result = (liste->size() == compareTo);
+  if (!result) {
+    DoLog(0) && (Log() << Verbose(0) << "Skipping because of the following atoms:" << std::endl);
+    for (LinkedCell::LinkedNodes::const_iterator iter = liste->begin();
+        iter != liste->end();
+        ++iter) {
+      DoLog(0) && (Log() << Verbose(0) << **iter << std::endl);
+    }
+  }
+  delete(liste);
+  return result;
+}
+
+/** Fills the empty space of the simulation box with water.
  * \param *filler molecule which the box is to be filled with
  * \param configuration contains box dimensions
  * \param distance[NDIM] distance between filling molecules in each direction
- * \param boundary length of boundary zone between molecule and filling mollecules
- * \param epsilon distance to surface which is not filled
+ * \param boundary length of boundary zone between molecule and filling molecules
  * \param RandAtomDisplacement maximum distance for random displacement per atom
  * \param RandMolDisplacement maximum distance for random displacement per filler molecule
- * \param DoRandomRotation true - do random rotiations, false - don't
- */
-void FillVoidWithMolecule(molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
+ * \param MinDistance minimum distance to boundary of domain and present molecules
+ * \param DoRandomRotation true - do random rotations, false - don't
+ */
+void FillVoidWithMolecule(
+    molecule *&filler,
+    config &configuration,
+    const double distance[NDIM],
+    const double boundary,
+    const double RandomAtomDisplacement,
+    const double RandomMolDisplacement,
+    const double MinDistance,
+    const bool DoRandomRotation
+    )
 {
   Info FunctionInfo(__func__);
@@ -949 +1052 @@
   RealSpaceMatrix Rotations;
   const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
-  Vector AtomTranslations;
   Vector FillerTranslations;
   Vector FillerDistance;
   Vector Inserter;
   double FillIt = false;
-  bond *Binder = NULL;
-  double phi[NDIM];
+  Vector firstInserter;
+  bool firstInsertion = true;
+  const Box &Domain = World::getInstance().getDomain();
   map<molecule *, LinkedCell *> LCList;
   std::vector<molecule *> List = World::getInstance().getAllMolecules();
@@ -968 +1071 @@
   // Center filler at origin
   filler->CenterEdge(&Inserter);
-  const int FillerCount = filler->getAtomCount();
+  //const int FillerCount = filler->getAtomCount();
   DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl);
   for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
@@ -996 +1099 @@
 
         // ... and rotation matrix
-        if (DoRandomRotation) {
-          for (int i=0;i<NDIM;i++) {
-            phi[i] = rand()/(RAND_MAX/(2.*M_PI));
-          }
-
-          Rotations.set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
-          Rotations.set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
-          Rotations.set(0,2,              cos(phi[1])*sin(phi[2])                                        );
-          Rotations.set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
-          Rotations.set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
-          Rotations.set(1,2,              sin(phi[1])                                                    );
-          Rotations.set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
-          Rotations.set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
-          Rotations.set(2,2,              cos(phi[1])*cos(phi[2])                                        );
+        if (DoRandomRotation)
+          Rotations.setRandomRotation();
+        else
+          Rotations.setIdentity();
+
+
+        // Check whether there is anything too close by and whether atom is outside of domain
+        FillIt = true;
+        for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {
+          FillIt = FillIt && isSpaceAroundPointVoid(
+              ListRunner->second,
+              (firstInsertion ? filler : NULL),
+              boundary,
+              CurrentPosition);
+          FillIt = FillIt && (Domain.isInside(CurrentPosition))
+              && ((Domain.DistanceToBoundary(CurrentPosition) - MinDistance) > -MYEPSILON);
+          if (!FillIt)
+            break;
         }
 
-        FillIt = true;
-        // go through all atoms
-        for(molecule::const_iterator iter = filler->begin(); iter !=filler->end();++iter){
-
-          // Check whether there is anything too close by
-          LinkedCell::LinkedNodes* liste = NULL;
-          for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {
-            liste = ListRunner->second->GetPointsInsideSphere(boundary, &CurrentPosition);
-            FillIt = FillIt && (liste->size() == 0);
-            delete(liste);
-            if (!FillIt)
-              break;
-          }
-        }
         // insert into Filling
         if (FillIt) {
           Inserter = CurrentPosition + FillerTranslations;
           DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is void point." << endl);
-          Filling = filler->CopyMolecule();
-          for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ++miter) {
-            if (DoRandomRotation) {
-              Vector temp = (*miter)->getPosition();
-              temp *= Rotations;
-              (*miter)->setPosition(temp);
-            }
-            // create atomic random translation vector ...
-            for (int i=0;i<NDIM;i++)
-              AtomTranslations[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
-            (*miter)->setPosition((*miter)->getPosition() + AtomTranslations);
+          // fill!
+          if (firstInsertion) { // use filler as first molecule
+            Filling = filler;
+            firstInsertion = false;
+            firstInserter = Inserter;
+          } else { // copy from filler molecule
+            Filling = filler->CopyMolecule();
+            RandomizeMoleculePositions(Filling, RandomAtomDisplacement, Rotations);
+            // translation
+            Filling->Translate(&Inserter);
+            // remove out-of-bounds atoms
+            RemoveAtomsOutsideDomain(Filling);
+            // TODO: Remove when World has no MoleculeListClass anymore
+            MolList->insert(Filling);
           }
-          Filling->Translate(&Inserter);
-          for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ) {
-            // check whether each atom is inside box
-            if (!World::getInstance().getDomain().isInside((*miter)->getPosition())) {
-              atom *Walker = *miter;
-              ++miter;
-              World::getInstance().destroyAtom(Walker);
-            } else {
-              ++miter;
-            }
-          }
-          MolList->insert(Filling);
         } else {
           DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is non-void point, within boundary or outside of MaxDistance." << endl);
@@ -1059 +1143 @@
         }
       }
-  // last one is replaced by the filler, as we need the original atoms contained therein!
-  filler->Translate(&Inserter);
-  MolList->erase(Filling);
-  for (molecule::iterator iter = Filling->begin(); !Filling->empty(); iter = Filling->begin()) {
-    atom *Walker = *iter;
-    Filling->erase(iter);
-    World::getInstance().destroyAtom(Walker);
-  }
-  World::getInstance().destroyMolecule(Filling);
+
+  // have we inserted any molecules?
+  if (firstInsertion) {
+    // If not remove filler
+    for(molecule::iterator miter = filler->begin(); !filler->empty(); miter = filler->begin()) {
+      atom *Walker = *miter;
+      filler->erase(Walker);
+      World::getInstance().destroyAtom(Walker);
+    }
+    World::getInstance().destroyMolecule(filler);
+  } else {
+    // otherwise translate and randomize to final position
+    if (DoRandomRotation)
+      Rotations.setRandomRotation();
+    else
+      Rotations.setIdentity();
+    RandomizeMoleculePositions(filler, RandomAtomDisplacement, Rotations);
+    // translation
+    filler->Translate(&firstInserter);
+    // remove out-of-bounds atoms
+    RemoveAtomsOutsideDomain(filler);
+  }
+
+  DoLog(0) && (Log() << Verbose(0) << MolList->ListOfMolecules.size() << " molecules have been inserted." << std::endl);
 
   for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) {