Index: src/Actions/MoleculeAction/FillVoidWithMoleculeAction.cpp
===================================================================
--- src/Actions/MoleculeAction/FillVoidWithMoleculeAction.cpp	(revision 0b15bbf617ce267e166ff4e7e1454ca420c1488f)
+++ src/Actions/MoleculeAction/FillVoidWithMoleculeAction.cpp	(revision 66fd499538327fa9d63d45ec481787c8650b9509)
@@ -28,4 +28,5 @@
 #include "World.hpp"
 
+#include "Descriptors/MoleculeIdDescriptor.hpp"
 #include "Parser/MpqcParser.hpp"
 #include "Parser/PcpParser.hpp"
@@ -35,4 +36,5 @@
 #include "Parser/FormatParserStorage.hpp"
 
+#include <algorithm>
 #include <iostream>
 #include <string>
@@ -47,11 +49,16 @@
 /** =========== define the function ====================== */
 Action::state_ptr MoleculeFillVoidWithMoleculeAction::performCall() {
-
   // obtain information
   getParametersfromValueStorage();
 
-  DoLog(1) && (Log() << Verbose(1) << "Filling Box with water molecules, lengths(" << params.lengths[0] << "," << params.lengths[1] << "," << params.lengths[2] << "), distances (" << params.distances[0] << "," << params.distances[1] << "," << params.distances[2] << "), DoRotate " << params.DoRotate << "." << endl);
+  if (!boost::filesystem::exists(params.fillername)) {
+    DoeLog(1) && (eLog() << Verbose(1) << "File with filler molecule " << params.fillername << " does not exist!" << endl);
+    return Action::failure;
+  }
+
+  DoLog(1) && (Log() << Verbose(1) << "Filling Box with water molecules, lengths(" << params.lengths[0] << "," << params.lengths[1] << "," << params.lengths[2] << "), distances (" << params.distances[0] << "," << params.distances[1] << "," << params.distances[2] << "), MinDistance " << params.MinDistance << ", DoRotate " << params.DoRotate << "." << endl);
   // construct water molecule
-  molecule *filler = World::getInstance().createMolecule();
+  std::vector<molecule *> presentmolecules = World::getInstance().getAllMolecules();
+//  DoLog(0) && (Log() << Verbose(0) << presentmolecules.size() << " molecules initially are present." << std::endl);
   std::string FilenameSuffix = params.fillername.string().substr(params.fillername.string().find_last_of('.')+1, params.fillername.string().length());
   ifstream input;
@@ -92,6 +99,10 @@
       break;
   }
-  World::MoleculeIterator iter = World::getInstance().getMoleculeIter();
-  for (; iter != World::getInstance().moleculeEnd(); ++iter)
+
+  // search the filler molecule that has been just parsed
+  molecule *filler = NULL;
+  for (World::MoleculeIterator iter = World::getInstance().getMoleculeIter();
+      iter != World::getInstance().moleculeEnd();
+      ++iter)
     filler = *iter; // get last molecule
   filler->SetNameFromFilename(params.fillername.string().c_str());
@@ -102,29 +113,64 @@
   for (int i=0;i<NDIM;i++)
     distance[i] = params.distances[i];
-  FillVoidWithMolecule(filler, *(World::getInstance().getConfig()), distance, params.lengths[0], params.lengths[1], params.lengths[2], params.DoRotate);
+  FillVoidWithMolecule(filler, *(World::getInstance().getConfig()), distance, params.lengths[0], params.lengths[1], params.lengths[2], params.MinDistance, params.DoRotate);
 
-  return Action::success;
+  // generate list of newly created molecules
+  // (we can in general remove more quickly from a list than a vector)
+  std::vector<molecule *> fillermolecules = World::getInstance().getAllMolecules();
+//  DoLog(0) && (Log() << Verbose(0) << fillermolecules.size() << " molecules are present." << std::endl);
+  std::list<molecule *> fillermolecules_list;
+  std::copy( fillermolecules.begin(),  fillermolecules.end(), std::back_inserter( fillermolecules_list ));
+//  DoLog(0) && (Log() << Verbose(0) << fillermolecules_list.size() << " molecules have been copied." << std::endl);
+  for (std::vector<molecule *>::const_iterator iter = presentmolecules.begin();
+      iter != presentmolecules.end();
+      ++iter) {
+    fillermolecules_list.remove(*iter);
+  }
+//  DoLog(0) && (Log() << Verbose(0) << fillermolecules_list.size() << " molecules left after removal." << std::endl);
+  fillermolecules.clear();
+  std::copy(fillermolecules_list.begin(), fillermolecules_list.end(), std::back_inserter( fillermolecules ));
+
+//  DoLog(0) && (Log() << Verbose(0) << fillermolecules.size() << " molecules have been inserted." << std::endl);
+
+  return Action::state_ptr(new MoleculeFillVoidWithMoleculeState(fillermolecules,params));
 }
 
 Action::state_ptr MoleculeFillVoidWithMoleculeAction::performUndo(Action::state_ptr _state) {
-//  MoleculeFillVoidWithMoleculeState *state = assert_cast<MoleculeFillVoidWithMoleculeState*>(_state.get());
+  MoleculeFillVoidWithMoleculeState *state = assert_cast<MoleculeFillVoidWithMoleculeState*>(_state.get());
 
-//  string newName = state->mol->getName();
-//  state->mol->setName(state->lastName);
+  MoleculeListClass *MolList = World::getInstance().getMolecules();
 
-  return Action::failure;
+  BOOST_FOREACH(molecule *_mol, state->fillermolecules) {
+    MolList->erase(_mol);
+    if ((_mol != NULL) && (!(World::getInstance().getAllMolecules(MoleculeById(_mol->getId()))).empty())) {
+      for(molecule::iterator iter = _mol->begin();
+          !_mol->empty();
+          iter = _mol->begin()) {
+        atom *Walker = *iter;
+        _mol->erase(iter);
+        World::getInstance().destroyAtom(Walker);
+      }
+      World::getInstance().destroyMolecule(_mol);
+    }
+  }
+
+  // as molecules and atoms from state are removed, we have to create a new one
+  std::vector<molecule *> fillermolecules;
+  return Action::state_ptr(new MoleculeFillVoidWithMoleculeState(fillermolecules,state->params));
 }
 
 Action::state_ptr MoleculeFillVoidWithMoleculeAction::performRedo(Action::state_ptr _state){
-  // Undo and redo have to do the same for this action
-  return performUndo(_state);
+  //MoleculeFillVoidWithMoleculeState *state = assert_cast<MoleculeFillVoidWithMoleculeState*>(_state.get());
+
+  return Action::failure;
+  //return Action::state_ptr(_state);
 }
 
 bool MoleculeFillVoidWithMoleculeAction::canUndo() {
-  return false;
+  return true;
 }
 
 bool MoleculeFillVoidWithMoleculeAction::shouldUndo() {
-  return false;
+  return true;
 }
 /** =========== end of function ====================== */
Index: src/Actions/MoleculeAction/FillVoidWithMoleculeAction.def
===================================================================
--- src/Actions/MoleculeAction/FillVoidWithMoleculeAction.def	(revision 0b15bbf617ce267e166ff4e7e1454ca420c1488f)
+++ src/Actions/MoleculeAction/FillVoidWithMoleculeAction.def	(revision 66fd499538327fa9d63d45ec481787c8650b9509)
@@ -13,12 +13,12 @@
 // ValueStorage by the token "Z" -> first column: int, Z, "Z"
 // "undefine" if no parameters are required, use (NODEFAULT) for each (undefined) default value
-#define paramtypes (boost::filesystem::path)(Vector)(Vector)(bool)
-#define paramtokens ("fill-void")("distances")("lengths")("DoRotate")
-#define paramdescriptions ("name of xyz file of filler molecule")("list of three of distances in space, one for each axis direction")("list of three of lengths in space, one for each axis direction")("whether to rotate or not")
-#undef paramdefaults
-#define paramreferences (fillername)(distances)(lengths)(DoRotate)
+#define paramtypes (boost::filesystem::path)(Vector)(Vector)(double)(bool)
+#define paramtokens ("fill-void")("distances")("lengths")("MinDistance")("DoRotate")
+#define paramdescriptions ("name of xyz file of filler molecule")("list of three of distances in space, one for each axis direction")("list of three of lengths in space, one for each axis direction")("minimum distance to boundary")("whether to rotate or not")
+#define paramdefaults (NODEFAULT)(NODEFAULT)(NODEFAULT)("0.")("0")
+#define paramreferences (fillername)(distances)(lengths)(MinDistance)(DoRotate)
 
-#undef statetypes
-#undef statereferences
+#define statetypes (std::vector<molecule *>)
+#define statereferences (fillermolecules)
 
 // some defines for all the names, you may use ACTION, STATE and PARAMS
Index: src/Box.cpp
===================================================================
--- src/Box.cpp	(revision 0b15bbf617ce267e166ff4e7e1454ca420c1488f)
+++ src/Box.cpp	(revision 66fd499538327fa9d63d45ec481787c8650b9509)
@@ -33,4 +33,6 @@
 
 #include "CodePatterns/Assert.hpp"
+#include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
 
 using namespace std;
@@ -220,9 +222,25 @@
 }
 
+double Box::DistanceToBoundary(const Vector &point) const
+{
+  std::map<double, Plane> DistanceSet;
+  std::vector<std::pair<Plane,Plane> > Boundaries = getBoundingPlanes();
+  for (int i=0;i<NDIM;++i) {
+    const double tempres1 = Boundaries[i].first.distance(point);
+    const double tempres2 = Boundaries[i].second.distance(point);
+    DistanceSet.insert( make_pair(tempres1, Boundaries[i].first) );
+    DoLog(1) && (Log() << Verbose(1) << "Inserting distance " << tempres1 << " and " << tempres2 << "." << std::endl);
+    DistanceSet.insert( make_pair(tempres2, Boundaries[i].second) );
+  }
+  ASSERT(!DistanceSet.empty(), "Box::DistanceToBoundary() - no distances in map!");
+  return (DistanceSet.begin())->first;
+}
+
 Shape Box::getShape() const{
   return transform(Cuboid(Vector(0,0,0),Vector(1,1,1)),(*M));
 }
 
-const Box::Conditions_t Box::getConditions(){
+const Box::Conditions_t Box::getConditions() const
+{
   return conditions;
 }
@@ -232,5 +250,6 @@
 }
 
-const vector<pair<Plane,Plane> >  Box::getBoundingPlanes(){
+const vector<pair<Plane,Plane> >  Box::getBoundingPlanes() const
+{
   vector<pair<Plane,Plane> > res;
   for(int i=0;i<NDIM;++i){
@@ -248,4 +267,5 @@
     res.push_back(make_pair(p1,p2));
   }
+  ASSERT(res.size() == 3, "Box::getBoundingPlanes() - does not have three plane pairs!");
   return res;
 }
Index: src/Box.hpp
===================================================================
--- src/Box.hpp	(revision 0b15bbf617ce267e166ff4e7e1454ca420c1488f)
+++ src/Box.hpp	(revision 66fd499538327fa9d63d45ec481787c8650b9509)
@@ -106,9 +106,15 @@
   double periodicDistance(const Vector &point1,const Vector &point2) const;
 
+  /**
+   * Calculates the minimum distance to the boundary of the periodic
+   * space defined by this box.
+   */
+  double DistanceToBoundary(const Vector &point) const;
+
   Shape getShape() const;
-  const Conditions_t getConditions();
+  const Conditions_t getConditions() const;
   void setCondition(int,BoundaryCondition_t);
 
-  const vector<pair<Plane,Plane> > getBoundingPlanes();
+  const vector<pair<Plane,Plane> > getBoundingPlanes() const;
 
   void setCuboid(const Vector&);
Index: src/LinearAlgebra/RealSpaceMatrix.cpp
===================================================================
--- src/LinearAlgebra/RealSpaceMatrix.cpp	(revision 0b15bbf617ce267e166ff4e7e1454ca420c1488f)
+++ src/LinearAlgebra/RealSpaceMatrix.cpp	(revision 66fd499538327fa9d63d45ec481787c8650b9509)
@@ -121,4 +121,24 @@
   set(2,2, cos(x)*cos(y));
 }
+
+void RealSpaceMatrix::setRandomRotation()
+{
+  double phi[NDIM];
+
+  for (int i=0;i<NDIM;i++) {
+    phi[i] = rand()/(RAND_MAX/(2.*M_PI));
+  }
+
+  set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
+  set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
+  set(0,2,              cos(phi[1])*sin(phi[2])                                        );
+  set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
+  set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
+  set(1,2,              sin(phi[1])                                                    );
+  set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
+  set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
+  set(2,2,              cos(phi[1])*cos(phi[2])                                        );
+}
+
 
 RealSpaceMatrix &RealSpaceMatrix::operator=(const RealSpaceMatrix &src)
Index: src/LinearAlgebra/RealSpaceMatrix.hpp
===================================================================
--- src/LinearAlgebra/RealSpaceMatrix.hpp	(revision 0b15bbf617ce267e166ff4e7e1454ca420c1488f)
+++ src/LinearAlgebra/RealSpaceMatrix.hpp	(revision 66fd499538327fa9d63d45ec481787c8650b9509)
@@ -72,4 +72,5 @@
    */
   void setRotation(const double x, const double y, const double z);
+  void setRandomRotation();
 
   /**
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 0b15bbf617ce267e166ff4e7e1454ca420c1488f)
+++ src/boundary.cpp	(revision 66fd499538327fa9d63d45ec481787c8650b9509)
@@ -927,17 +927,120 @@
 };
 
-/** 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,11 +1052,11 @@
   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,5 +1071,5 @@
   // 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,62 +1099,43 @@
 
         // ... 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,13 +1143,28 @@
         }
       }
-  // 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()) {
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 0b15bbf617ce267e166ff4e7e1454ca420c1488f)
+++ src/boundary.hpp	(revision 66fd499538327fa9d63d45ec481787c8650b9509)
@@ -41,5 +41,5 @@
 double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
 molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
-void FillVoidWithMolecule(molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
+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);
     void FindConvexBorder(const molecule* const mol, Boundaries *BoundaryPts, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename);
 Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol);
