Index: src/Actions/Values.cpp
===================================================================
--- src/Actions/Values.cpp	(revision c96423eec2163f37cc351c88bfd115e0a0446366)
+++ src/Actions/Values.cpp	(revision 2a0271161126f3b47c3cd19c934d5f9e68593ecf)
@@ -41,5 +41,5 @@
   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)+".");
 
Index: src/Box.cpp
===================================================================
--- src/Box.cpp	(revision c96423eec2163f37cc351c88bfd115e0a0446366)
+++ src/Box.cpp	(revision 2a0271161126f3b47c3cd19c934d5f9e68593ecf)
@@ -169,4 +169,18 @@
     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)));
Index: src/Box.hpp
===================================================================
--- src/Box.hpp	(revision c96423eec2163f37cc351c88bfd115e0a0446366)
+++ src/Box.hpp	(revision 2a0271161126f3b47c3cd19c934d5f9e68593ecf)
@@ -83,7 +83,13 @@
 
   /**
-   * 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,4 +152,5 @@
   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;
Index: src/LinkedCell/LinkedCell_Model.cpp
===================================================================
--- src/LinkedCell/LinkedCell_Model.cpp	(revision c96423eec2163f37cc351c88bfd115e0a0446366)
+++ src/LinkedCell/LinkedCell_Model.cpp	(revision 2a0271161126f3b47c3cd19c934d5f9e68593ecf)
@@ -288,4 +288,7 @@
 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));
Index: src/LinkedCell/unittests/stubs/ObserverBoxStub.cpp
===================================================================
--- src/LinkedCell/unittests/stubs/ObserverBoxStub.cpp	(revision c96423eec2163f37cc351c88bfd115e0a0446366)
+++ src/LinkedCell/unittests/stubs/ObserverBoxStub.cpp	(revision 2a0271161126f3b47c3cd19c934d5f9e68593ecf)
@@ -70,4 +70,9 @@
 }
 
+bool Box::isValid(const Vector &p) const
+{
+  return true;
+}
+
 void Box::setM(RealSpaceMatrix _M)
 {
Index: src/Tesselation/boundary.cpp
===================================================================
--- src/Tesselation/boundary.cpp	(revision c96423eec2163f37cc351c88bfd115e0a0446366)
+++ src/Tesselation/boundary.cpp	(revision 2a0271161126f3b47c3cd19c934d5f9e68593ecf)
@@ -1175,5 +1175,5 @@
               boundary,
               CurrentPosition);
-          FillIt = FillIt && (Domain.isInside(CurrentPosition))
+          FillIt = FillIt && (Domain.isValid(CurrentPosition))
               && ((Domain.DistanceToBoundary(CurrentPosition) - MinDistance) > -MYEPSILON);
           if (!FillIt)
@@ -1240,4 +1240,217 @@
 };
 
+
+/** 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
Index: src/UIElements/CommandLineUI/Query/VectorCommandLineQuery.cpp
===================================================================
--- src/UIElements/CommandLineUI/Query/VectorCommandLineQuery.cpp	(revision c96423eec2163f37cc351c88bfd115e0a0446366)
+++ src/UIElements/CommandLineUI/Query/VectorCommandLineQuery.cpp	(revision 2a0271161126f3b47c3cd19c934d5f9e68593ecf)
@@ -41,5 +41,5 @@
     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;
Index: src/UIElements/CommandLineUI/Query/VectorsCommandLineQuery.cpp
===================================================================
--- src/UIElements/CommandLineUI/Query/VectorsCommandLineQuery.cpp	(revision c96423eec2163f37cc351c88bfd115e0a0446366)
+++ src/UIElements/CommandLineUI/Query/VectorsCommandLineQuery.cpp	(revision 2a0271161126f3b47c3cd19c934d5f9e68593ecf)
@@ -42,5 +42,5 @@
     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
Index: src/UIElements/TextUI/Query/VectorTextQuery.cpp
===================================================================
--- src/UIElements/TextUI/Query/VectorTextQuery.cpp	(revision c96423eec2163f37cc351c88bfd115e0a0446366)
+++ src/UIElements/TextUI/Query/VectorTextQuery.cpp	(revision 2a0271161126f3b47c3cd19c934d5f9e68593ecf)
@@ -67,5 +67,5 @@
 
   // check vector
-  return World::getInstance().getDomain().isInside(tmp);
+  return World::getInstance().getDomain().isValid(tmp);
 }
 
Index: src/UIElements/TextUI/Query/VectorsTextQuery.cpp
===================================================================
--- src/UIElements/TextUI/Query/VectorsTextQuery.cpp	(revision c96423eec2163f37cc351c88bfd115e0a0446366)
+++ src/UIElements/TextUI/Query/VectorsTextQuery.cpp	(revision 2a0271161126f3b47c3cd19c934d5f9e68593ecf)
@@ -71,5 +71,5 @@
         temp[counter++] = coord;
       }
-      if (World::getInstance().getDomain().isInside(temp))
+      if (World::getInstance().getDomain().isValid(temp))
         tmp.push_back(temp);
       olditerspace = vectoriter;
