Changeset 2a0271 for src

Timestamp:

Feb 22, 2012, 11:29:22 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:

7e3065

Parents:

c96423

git-author:

Frederik Heber <heber@…> (01/18/12 18:13:52)

git-committer:

Frederik Heber <heber@…> (02/22/12 11:29:22)

Message:

New function Box::isValid() to be discerned from ::isInside().

• isValid checks with and isInside without boundary conditions.
• added Box::isValid() check to LinkedCell_Model::getIndexToVector().
• TESTFIX: regression test Atoms/Add needed Ignore BCs.

Location:

src

Files:

• Actions/Values.cpp (modified) (1 diff)
• Box.cpp (modified) (1 diff)
• Box.hpp (modified) (2 diffs)
• LinkedCell/LinkedCell_Model.cpp (modified) (1 diff)
• LinkedCell/unittests/stubs/ObserverBoxStub.cpp (modified) (1 diff)
• Tesselation/boundary.cpp (modified) (2 diffs)
• UIElements/CommandLineUI/Query/VectorCommandLineQuery.cpp (modified) (1 diff)
• UIElements/CommandLineUI/Query/VectorsCommandLineQuery.cpp (modified) (1 diff)
• UIElements/TextUI/Query/VectorTextQuery.cpp (modified) (1 diff)
• UIElements/TextUI/Query/VectorsTextQuery.cpp (modified) (1 diff)

src/Actions/Values.cpp

@@ -41 +41 @@
   static_cast<Vector>(returnVector) = Vector(vector); // under its hood it's still a Vector
 
-  ASSERT(_box.isInside(returnVector),
+  ASSERT(_box.isValid(returnVector),
       "VectorValue::toBoxVector() - vector "+toString(toVector())+" does not lie within box "+toString(_box)+".");
 

src/Box.cpp

@@ -169 +169 @@
     result = result &&
               ((conditions[i] == BoundaryConditions::Ignore) ||
+               ((tester[i] >= -MYEPSILON) &&
+               ((tester[i] - 1.) < MYEPSILON)));
+
+  return result;
+}
+
+bool Box::isValid(const Vector &point) const
+{
+  bool result = true;
+  Vector tester = translateOut(point);
+
+  for(int i=0;i<NDIM;i++)
+    result = result &&
+              ((conditions[i] != BoundaryConditions::Ignore) ||
                ((tester[i] >= -MYEPSILON) &&
                ((tester[i] - 1.) < MYEPSILON)));

src/Box.hpp

@@ -83 +83 @@
 
   /**
-   * Checks whether a given vector is inside the box.
+   * Checks whether a given vector is inside the box not accounting for boundary conditions.
    */
   bool isInside(const Vector &point) const;
+
+  /**
+   * Checks whether a given vector is inside the box under the given boundary conditions.
+   */
+  bool isValid(const Vector &point) const;
+
 
   /**
@@ -146 +152 @@
   void internal_explode(const Vector &point,int n) const;
 
+
   //!> Internal vector list for exploding vectors used in Box::internal_explode().
   mutable VECTORSET(std::vector) internal_list;

src/LinkedCell/LinkedCell_Model.cpp

@@ -288 +288 @@
 const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getIndexToVector(const Vector &position) const
 {
+  ASSERT(domain.isValid(position),
+      "LinkedCell::LinkedCell_Model::getIndexToVector() - specified position "+toString(position)
+      +"is not valid in the current domain.");
   tripleIndex index;
   Vector x(Partition*domain.enforceBoundaryConditions(position));

src/LinkedCell/unittests/stubs/ObserverBoxStub.cpp

@@ -70 +70 @@
 }
 
+bool Box::isValid(const Vector &p) const
+{
+  return true;
+}
+
 void Box::setM(RealSpaceMatrix _M)
 {

src/Tesselation/boundary.cpp

@@ -1175 +1175 @@
               boundary,
               CurrentPosition);
-          FillIt = FillIt && (Domain.isInside(CurrentPosition))
+          FillIt = FillIt && (Domain.isValid(CurrentPosition))
               && ((Domain.DistanceToBoundary(CurrentPosition) - MinDistance) > -MYEPSILON);
           if (!FillIt)
@@ -1240 +1240 @@
 };
 
+
+/** Fills the empty space around other molecules' surface of the simulation box with a filler.
+ *
+ * Note that we use \a FindNonConvexBorder to determine the surface of the found molecules.
+ * There, we use a radius of twice the given \a boundary.
+ *
+ * \param *out output stream for debugging
+ * \param *List list of molecules already present in box
+ * \param *TesselStruct contains tesselated surface
+ * \param *filler molecule which the box is to be filled with
+ * \param configuration contains box dimensions
+ * \param MaxSurfaceDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
+ * \param distance[NDIM] distance between filling molecules in each direction
+ * \param boundary length of boundary zone between molecule and filling molecules
+ * \param MinDistance distance to boundary of domain which is not filled
+ * \param RandAtomDisplacement maximum distance for random displacement per atom
+ * \param RandMolDisplacement maximum distance for random displacement per filler molecule
+ * \param DoRandomRotation true - do random rotations, false - don't
+ */
+void FillBoxWithMolecule(
+    MoleculeListClass *MolList,
+    molecule *filler,
+    config &configuration,
+    const double MaxSurfaceDistance,
+    const double distance[NDIM],
+    const double boundary,
+    const double MinDistance,
+    const double RandomAtomDisplacement,
+    const double RandomMolDisplacement,
+    const bool DoRandomRotation)
+{
+  Info FunctionInfo(__func__);
+  molecule *Filling = NULL;
+  Vector CurrentPosition;
+  int N[NDIM];
+  int n[NDIM];
+  const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
+  RealSpaceMatrix Rotations;
+  const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
+  Vector FillerTranslations;
+  Vector FillerDistance;
+  Vector Inserter;
+  double FillIt = false;
+  Vector firstInserter;
+  bool firstInsertion = true;
+  const Box &Domain = World::getInstance().getDomain();
+  map<molecule *, LinkedCell_deprecated *> LCList;
+  std::vector<molecule *> List = World::getInstance().getAllMolecules();
+  map<molecule *, Tesselation *> TesselStruct;
+
+  for (MoleculeList::iterator ListRunner = MolList->ListOfMolecules.begin(); ListRunner != MolList->ListOfMolecules.end(); ListRunner++)
+    if ((*ListRunner)->getAtomCount() > 0) {
+      LOG(1, "Pre-creating linked cell lists for molecule " << *ListRunner << ".");
+      PointCloudAdaptor< molecule > cloud(*ListRunner, (*ListRunner)->name);
+      LCList[(*ListRunner)] = new LinkedCell_deprecated(cloud, 4.*boundary); // get linked cell list
+      LOG(1, "Pre-creating tesselation for molecule " << *ListRunner << ".");
+      TesselStruct[(*ListRunner)] = NULL;
+      FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell_deprecated *&)LCList[(*ListRunner)], 2.*boundary, NULL);
+    }
+
+  // Center filler at origin
+  filler->CenterEdge(&Inserter);
+  const int FillerCount = filler->getAtomCount();
+  LOG(2, "INFO: Filler molecule has the following bonds:");
+  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
+    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
+    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
+        BondRunner != ListOfBonds.end();
+        ++BondRunner) {
+      if ((*BondRunner)->leftatom == *AtomRunner)
+        LOG(2, "  " << *(*BondRunner));
+    }
+  }
+
+  atom * CopyAtoms[FillerCount];
+
+  setVerbosity(4);
+
+  // calculate filler grid in [0,1]^3
+  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
+  for(int i=0;i<NDIM;i++)
+    N[i] = (int) ceil(1./FillerDistance[i]);
+  LOG(1, "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << ".");
+
+  // initialize seed of random number generator to current time
+  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
+  const double rng_min = random.min();
+  const double rng_max = random.max();
+  //srand ( time(NULL) );
+
+  // go over [0,1]^3 filler grid
+  for (n[0] = 0; n[0] < N[0]; n[0]++)
+    for (n[1] = 0; n[1] < N[1]; n[1]++)
+      for (n[2] = 0; n[2] < N[2]; n[2]++) {
+        // calculate position of current grid vector in untransformed box
+        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
+        // create molecule random translation vector ...
+        for (int i=0;i<NDIM;i++)
+          FillerTranslations[i] = RandomMolDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
+        LOG(1, "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << ".");
+
+        // ... and rotation matrix
+        if (DoRandomRotation)
+          setRandomRotation(Rotations);
+        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_deprecated *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {
+          // check whether its void
+          FillIt = FillIt && isSpaceAroundPointVoid(
+              ListRunner->second,
+              (firstInsertion ? filler : NULL),
+              boundary,
+              CurrentPosition);
+          if (!FillIt) {
+            LOG(2, "REJECT: Position at " << Inserter << " is non-void.");
+            break;
+          }
+          // check whether inside domain
+          FillIt = FillIt && (Domain.isValid(CurrentPosition));
+          if (!FillIt) {
+            LOG(2, "REJECT: Position at " << CurrentPosition << " is "
+                << distance << ", hence outside of domain.");
+            break;
+          }
+          // check minimum distance to boundary
+          const double distance = (Domain.DistanceToBoundary(CurrentPosition) - MinDistance);
+          FillIt = FillIt && (distance > -MYEPSILON);
+          if (!FillIt) {
+            LOG(2, "REJECT: Position at " << CurrentPosition << " is " << distance << ", less than "
+                << MinDistance << " hence, too close to boundary.");
+            break;
+          }
+        }
+        // Check whether point is in- or outside of tesselations
+        if (FillIt) {
+          for (MoleculeList::iterator ListRunner = MolList->ListOfMolecules.begin(); ListRunner != MolList->ListOfMolecules.end(); ListRunner++) {
+            // get linked cell list
+            if (TesselStruct[(*ListRunner)] != NULL) {
+              const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
+              LOG(2, "INFO: Distance to surface is " << distance << ".");
+              FillIt = FillIt && ((distance == -1.) || (distance > boundary)) && ((MaxSurfaceDistance < 0) || (MaxSurfaceDistance > distance));
+              if (!FillIt) {
+                LOG(2, "REJECT: Position at " << CurrentPosition << " is in distance of " << distance
+                    << " to a surface, less than " << MaxSurfaceDistance  << " hence, too close.");
+                break;
+              }
+            }
+          }
+        }
+
+        // insert into Filling
+        if (FillIt) {
+          Inserter = CurrentPosition + FillerTranslations;
+          LOG(2, "ACCEPT: Position at " << CurrentPosition << " is void point.");
+          // fill!
+          Filling = filler->CopyMolecule();
+          RandomizeMoleculePositions(Filling, RandomAtomDisplacement, Rotations, random);
+          // translation
+          Filling->Translate(&Inserter);
+          // remove out-of-bounds atoms
+          const bool status = RemoveAtomsOutsideDomain(Filling);
+          if ((firstInsertion) && (!status)) { // no atom has been removed
+            // remove copied atoms and molecule again
+            Filling->removeAtomsinMolecule();
+            World::getInstance().destroyMolecule(Filling);
+            // and mark is final filler position
+            Filling = filler;
+            firstInsertion = false;
+            firstInserter = Inserter;
+          } else {
+            // TODO: Remove when World has no MoleculeListClass anymore
+            if (Filling)
+              MolList->insert(Filling);
+          }
+        } else {
+          LOG(2, "REJECT: Position at " << CurrentPosition << " is non-void point, within boundary or outside of MaxSurfaceDistance.");
+          continue;
+        }
+      }
+
+  // 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;
+      World::getInstance().destroyAtom(Walker);
+    }
+    World::getInstance().destroyMolecule(filler);
+  } else {
+    // otherwise translate and randomize to final position
+    if (DoRandomRotation)
+      setRandomRotation(Rotations);
+    else
+      Rotations.setIdentity();
+    RandomizeMoleculePositions(filler, RandomAtomDisplacement, Rotations, random);
+    // translation
+    filler->Translate(&firstInserter);
+    // remove out-of-bounds atoms
+    RemoveAtomsOutsideDomain(filler);
+  }
+
+  LOG(0, MolList->ListOfMolecules.size() << " molecules have been inserted.");
+
+  for (std::map<molecule *, LinkedCell_deprecated *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) {
+    delete ListRunner->second;
+    LCList.erase(ListRunner);
+  }
+};
+
 /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
  * \param *out output stream for debugging

src/UIElements/CommandLineUI/Query/VectorCommandLineQuery.cpp

@@ -41 +41 @@
     temp = CommandLineParser::getInstance().vm[getTitle()].as< VectorValue >();
     tmp = temp.toVector();
-    if ((check) && (!World::getInstance().getDomain().isInside(tmp))) {
+    if ((check) && (!World::getInstance().getDomain().isValid(tmp))) {
       ELOG(1, "Vector " << tmp << " would be outside of box domain.");
       return false;

src/UIElements/CommandLineUI/Query/VectorsCommandLineQuery.cpp

@@ -42 +42 @@
     for(std::vector<VectorValue>::iterator iter = temporary.begin(); iter != temporary.end(); ++iter) {
       temp = (*iter).toVector();
-      if ((!check) || (World::getInstance().getDomain().isInside(temp)))
+      if ((!check) || (World::getInstance().getDomain().isValid(temp)))
         tmp.push_back(temp);
       else

src/UIElements/TextUI/Query/VectorTextQuery.cpp

@@ -67 +67 @@
 
   // check vector
-  return World::getInstance().getDomain().isInside(tmp);
+  return World::getInstance().getDomain().isValid(tmp);
 }
 

src/UIElements/TextUI/Query/VectorsTextQuery.cpp

@@ -71 +71 @@
         temp[counter++] = coord;
       }
-      if (World::getInstance().getDomain().isInside(temp))
+      if (World::getInstance().getDomain().isValid(temp))
         tmp.push_back(temp);
       olditerspace = vectoriter;