Index: src/Actions/FillAction/SuspendInMoleculeAction.cpp
===================================================================
--- src/Actions/FillAction/SuspendInMoleculeAction.cpp	(revision 55736be799388502de9a57ac542d215e5ce0565a)
+++ src/Actions/FillAction/SuspendInMoleculeAction.cpp	(revision 55736be799388502de9a57ac542d215e5ce0565a)
@@ -0,0 +1,333 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2014 Frederik Heber. All rights reserved.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * SuspendInMoleculeAction.cpp
+ *
+ *  Created on: Sep 05, 2014
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "Actions/UndoRedoHelpers.hpp"
+#include "Atom/atom.hpp"
+#include "Atom/AtomicInfo.hpp"
+#include "Atom/CopyAtoms/CopyAtoms_withBonds.hpp"
+#include "Bond/BondInfo.hpp"
+#include "CodePatterns/Log.hpp"
+#include "Descriptors/MoleculeOrderDescriptor.hpp"
+#include "Element/element.hpp"
+#include "Filling/Cluster.hpp"
+#include "Filling/Filler.hpp"
+#include "Filling/Preparators/BoxFillerPreparator.hpp"
+#include "LinkedCell/linkedcell.hpp"
+#include "LinkedCell/PointCloudAdaptor.hpp"
+#include "molecule.hpp"
+#include "MoleculeListClass.hpp"
+#include "Parser/FormatParserInterface.hpp"
+#include "Parser/FormatParserStorage.hpp"
+#include "Tesselation/boundary.hpp"
+#include "Tesselation/tesselation.hpp"
+#include "World.hpp"
+
+#include <algorithm>
+#include<gsl/gsl_poly.h>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "Actions/FillAction/SuspendInMoleculeAction.hpp"
+
+using namespace MoleCuilder;
+
+static double calculateMass(const molecule &_mol)
+{
+  // sum up the atomic masses
+  const double mass = _mol.getAtomSet().totalMass();
+  LOG(2, "DEBUG: Molecule "+_mol.getName()+"'s summed mass is "
+      << setprecision(10) << mass << " atomicmassunit.");
+  return mass;
+}
+
+static double calculateEnvelopeVolume(
+    molecule &_mol,
+    std::vector<double> &_diameters)
+{
+  const bool IsAngstroem = true;
+  class Tesselation *TesselStruct = NULL;
+
+  Boundaries *BoundaryPoints = GetBoundaryPoints(&_mol, TesselStruct);
+  const double * diameters =
+      GetDiametersOfCluster(BoundaryPoints, &_mol, TesselStruct, IsAngstroem);
+  std::copy(&diameters[0], &diameters[3], _diameters.begin());
+  delete diameters;
+  PointCloudAdaptor< molecule > cloud(&_mol, _mol.getName());
+  LinkedCell_deprecated *LCList = new LinkedCell_deprecated(cloud, 10.);
+  FindConvexBorder(&_mol, BoundaryPoints, TesselStruct, (const LinkedCell_deprecated *&)LCList, NULL);
+  delete (LCList);
+  delete[] BoundaryPoints;
+
+  // some preparations beforehand
+  const double volume = TesselStruct->getVolumeOfConvexEnvelope(IsAngstroem);
+
+  delete TesselStruct;
+
+  LOG(2, "DEBUG: Molecule "+_mol.getName()+"'s volume is "
+      << setprecision(10) << volume << " angstrom^3.");
+
+  return volume;
+}
+
+// and construct the stuff
+#include "SuspendInMoleculeAction.def"
+#include "Action_impl_pre.hpp"
+/** =========== define the function ====================== */
+ActionState::ptr FillSuspendInMoleculeAction::performCall() {
+  typedef std::vector<atom*> AtomVector;
+
+  // get the filler molecule
+  std::vector<AtomicInfo> movedatoms;
+  molecule *filler = NULL;
+  {
+    const std::vector< molecule *> molecules = World::getInstance().getSelectedMolecules();
+    if (molecules.size() != 1) {
+      STATUS("No exactly one molecule selected, aborting,");
+      return Action::failure;
+    }
+    filler = *(molecules.begin());
+  }
+  for(molecule::const_iterator iter = filler->begin(); iter != filler->end(); ++iter)
+    movedatoms.push_back( AtomicInfo(*(*iter)) );
+  LOG(1, "INFO: Chosen molecule has " << filler->size() << " atoms.");
+
+  // center filler's tip at origin
+  filler->CenterEdge();
+
+  /// first we need to calculate some volumes and masses
+  std::vector<molecule *> molecules = World::getInstance().getAllMolecules();
+  double totalmass = 0.;
+  const bool IsAngstroem = true;
+  Vector BoxLengths;
+  double clustervolume = 0.;
+  std::vector<double> GreatestDiameter(NDIM, 0.);
+  for (std::vector<molecule *>::iterator iter = molecules.begin();
+      iter != molecules.end(); ++iter)
+  {
+    molecule & mol = **iter;
+    const double mass = calculateMass(mol);
+    totalmass += mass;
+    std::vector<double> diameters(NDIM, 0.);
+    const double volume = calculateEnvelopeVolume(mol, diameters);
+    clustervolume += volume;
+    for (size_t i=0;i<NDIM;++i)
+      GreatestDiameter[i] = std::max(GreatestDiameter[i], diameters[i]);
+  }
+  LOG(1, "INFO: The summed mass is " << setprecision(10)
+      << totalmass << " atomicmassunit.");
+  LOG(1, "INFO: The average density is " << setprecision(10)
+      << totalmass / clustervolume << " atomicmassunit/"
+      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+
+  // calculate maximum solvent density
+  std::vector<double> fillerdiameters(NDIM, 0.);
+  const double solventdensity =
+      calculateMass(*filler) / calculateEnvelopeVolume(*filler, fillerdiameters);
+
+  std::vector<unsigned int> counts(3, 0);
+  Vector offset(.5,.5,.5);
+
+  /// solve cubic polynomial
+  double cellvolume = 0.;
+  LOG(1, "Solving equidistant suspension in water problem ...");
+  cellvolume = (totalmass / solventdensity
+      - (totalmass / clustervolume)) / (params.density.get() - 1.);
+  LOG(1, "Cellvolume needed for a density of " << params.density.get()
+      << " g/cm^3 is " << cellvolume << " angstroem^3.");
+
+  const double minimumvolume =
+      (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
+  LOG(1, "Minimum volume of the convex envelope contained in a rectangular box is "
+      << minimumvolume << "angstrom^3.");
+
+  if (minimumvolume > cellvolume) {
+    ELOG(1, "The containing box already has a greater volume than the envisaged cell volume!");
+    LOG(0, "Setting Box dimensions to minimum possible, the greatest diameters.");
+    for (int i = 0; i < NDIM; i++)
+      BoxLengths[i] = GreatestDiameter[i];
+//    mol->CenterEdge();
+  } else {
+    BoxLengths[0] = GreatestDiameter[0] + GreatestDiameter[1] + GreatestDiameter[2];
+    BoxLengths[1] = GreatestDiameter[0] * GreatestDiameter[1]
+        + GreatestDiameter[0] * GreatestDiameter[2]
+        + GreatestDiameter[1] * GreatestDiameter[2];
+    BoxLengths[2] = minimumvolume - cellvolume;
+    double x0 = 0.;
+    double x1 = 0.;
+    double x2 = 0.;
+    // for cubic polynomial there are either 1 or 3 unique solutions
+    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x0, &x1, &x2) == 1)
+      LOG(0, "RESULT: The resulting spacing is: " << x0 << " .");
+    else {
+      LOG(0, "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " .");
+      std::swap(x0, x2); // sorted in ascending order
+    }
+
+    cellvolume = 1.;
+    for (size_t i = 0; i < NDIM; ++i) {
+      BoxLengths[i] = x0 + GreatestDiameter[i];
+      cellvolume *= BoxLengths[i];
+    }
+
+    // set new box dimensions
+    LOG(0, "Translating to box with these boundaries.");
+    {
+      RealSpaceMatrix domain;
+      for(size_t i =0; i<NDIM;++i)
+        domain.at(i,i) = BoxLengths[i];
+      World::getInstance().setDomain(domain);
+    }
+//    mol->CenterInBox();
+  }
+
+  // update Box of atoms by boundary
+  {
+    RealSpaceMatrix domain;
+    for(size_t i =0; i<NDIM;++i)
+      domain.at(i,i) = BoxLengths[i];
+    World::getInstance().setDomain(domain);
+  }
+  LOG(0, "RESULT: The resulting cell dimensions are: " << BoxLengths[0] << " and " << BoxLengths[1] << " and " << BoxLengths[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+
+  // prepare the filler preparator
+  BoxFillerPreparator filler_preparator(filler);
+  filler_preparator.addVoidPredicate(params.mindistance.get());
+  filler_preparator.addRandomInserter(
+      params.RandAtomDisplacement.get(),
+      params.RandMoleculeDisplacement.get(),
+      params.DoRotate.get());
+  filler_preparator.addCubeMesh(
+      counts,
+      offset,
+      World::getInstance().getDomain().getM());
+  if (!filler_preparator()) {
+    STATUS("Filler was not fully constructed.");
+    return Action::failure;
+  }
+
+  // use filler
+  bool successflag = false;
+  FillSuspendInMoleculeState *UndoState = NULL;
+  {
+    // fill
+    Filler *fillerFunction = filler_preparator.obtainFiller();
+    // TODO: When molecule::getBoundingSphere() does not use a sphere anymore,
+    // we need to check whether we rotate the molecule randomly. For this to
+    // work we need a sphere!
+    const Shape s = filler->getBoundingSphere(params.RandAtomDisplacement.get());
+    ClusterInterface::Cluster_impl cluster( new Cluster(filler->getAtomIds(), s) );
+    CopyAtoms_withBonds copyMethod;
+    Filler::ClusterVector_t ClonedClusters;
+    successflag = (*fillerFunction)(copyMethod, cluster, ClonedClusters);
+    delete fillerFunction;
+
+    // append each cluster's atoms to clonedatoms (however not selected ones)
+    std::vector<const atom *> clonedatoms;
+    std::vector<AtomicInfo> clonedatominfos;
+    for (Filler::ClusterVector_t::const_iterator iter = ClonedClusters.begin();
+        iter != ClonedClusters.end(); ++iter) {
+      const AtomIdSet &atoms = (*iter)->getAtomIds();
+      clonedatoms.reserve(clonedatoms.size()+atoms.size());
+      for (AtomIdSet::const_iterator atomiter = atoms.begin(); atomiter != atoms.end(); ++atomiter)
+        if (!filler->containsAtom(*atomiter)) {
+          clonedatoms.push_back( *atomiter );
+          clonedatominfos.push_back( AtomicInfo(*(*atomiter)) );
+        }
+    }
+    std::vector< BondInfo > clonedbonds;
+    StoreBondInformationFromAtoms(clonedatoms, clonedbonds);
+    LOG(2, "DEBUG: There are " << clonedatominfos.size() << " newly created atoms.");
+
+    if (!successflag) {
+      STATUS("Insertion failed, removing inserted clusters, translating original one back");
+      RemoveAtomsFromAtomicInfo(clonedatominfos);
+      clonedatoms.clear();
+      SetAtomsFromAtomicInfo(movedatoms);
+    } else {
+      std::vector<Vector> MovedToVector(filler->size(), zeroVec);
+      std::transform(filler->begin(), filler->end(), MovedToVector.begin(),
+          boost::bind(&AtomInfo::getPosition, _1) );
+      UndoState = new FillSuspendInMoleculeState(clonedatominfos,clonedbonds,movedatoms,MovedToVector,params);
+    }
+  }
+
+  if (successflag)
+    return ActionState::ptr(UndoState);
+  else {
+    return Action::failure;
+  }
+}
+
+ActionState::ptr FillSuspendInMoleculeAction::performUndo(ActionState::ptr _state) {
+  FillSuspendInMoleculeState *state = assert_cast<FillSuspendInMoleculeState*>(_state.get());
+
+  // remove all created atoms
+  RemoveAtomsFromAtomicInfo(state->clonedatoms);
+  // add the original cluster
+  SetAtomsFromAtomicInfo(state->movedatoms);
+
+  return ActionState::ptr(_state);
+}
+
+ActionState::ptr FillSuspendInMoleculeAction::performRedo(ActionState::ptr _state){
+  FillSuspendInMoleculeState *state = assert_cast<FillSuspendInMoleculeState*>(_state.get());
+
+  // place filler cluster again at new spot
+  ResetAtomPosition(state->movedatoms, state->MovedToVector);
+
+  // re-create all clusters
+  bool statusflag = AddAtomsFromAtomicInfo(state->clonedatoms);
+
+  // re-create the bonds
+  if (statusflag)
+    AddBondsFromBondInfo(state->clonedbonds);
+  if (statusflag)
+    return ActionState::ptr(_state);
+  else {
+    STATUS("Failed re-adding filled in atoms.");
+    return Action::failure;
+  }
+}
+
+bool FillSuspendInMoleculeAction::canUndo() {
+  return false;
+}
+
+bool FillSuspendInMoleculeAction::shouldUndo() {
+  return false;
+}
+/** =========== end of function ====================== */
Index: src/Actions/FillAction/SuspendInMoleculeAction.def
===================================================================
--- src/Actions/FillAction/SuspendInMoleculeAction.def	(revision 55736be799388502de9a57ac542d215e5ce0565a)
+++ src/Actions/FillAction/SuspendInMoleculeAction.def	(revision 55736be799388502de9a57ac542d215e5ce0565a)
@@ -0,0 +1,52 @@
+/*
+ * SuspendInMoleculeAction.def
+ *
+ *  Created on: Aug 26, 2010
+ *      Author: heber
+ */
+
+// all includes and forward declarations necessary for non-integral types below
+#include <limits>
+#include <vector>
+#include "Atom/AtomicInfo.hpp"
+#include "Bond/BondInfo.hpp"
+#include "types.hpp"
+
+#include "LinearAlgebra/defs.hpp"
+#include "Parameters/Validators/DummyValidator.hpp"
+#include "Parameters/Validators/RangeValidator.hpp"
+#include "Parameters/Validators/STLVectorValidator.hpp"
+#include "Parameters/Validators/Specific/BoxLengthValidator.hpp"
+#include "Parameters/Validators/Specific/VectorZeroOneComponentsValidator.hpp"
+
+// i.e. there is an integer with variable name Z that can be found in
+// ValueStorage by the token "Z" -> first column: int, Z, "Z"
+// "undefine" if no parameters are required, use (NOPARAM_DEFAULT) for each (undefined) default value
+#define paramtypes (double)(double)(double)(double)(bool)
+#define paramtokens ("density")("min-distance")("random-atom-displacement")("random-molecule-displacement")("DoRotate")
+#define paramdescriptions ("desired density for the total domain")("minimum distance of water molecules to present atoms")("magnitude of random atom displacement")("magnitude of random molecule displacement")("whether to rotate or not")
+#define paramdefaults (PARAM_DEFAULT(1.))(PARAM_DEFAULT(1.))(PARAM_DEFAULT(0.))(PARAM_DEFAULT(0.))(PARAM_DEFAULT(false))
+#define paramreferences (density)(mindistance)(RandAtomDisplacement)(RandMoleculeDisplacement)(DoRotate)
+#define paramvalids \
+(RangeValidator< double >(0., std::numeric_limits<double>::max())) \
+(BoxLengthValidator()) \
+(BoxLengthValidator()) \
+(BoxLengthValidator()) \
+(DummyValidator< bool >())
+
+#define statetypes (std::vector<AtomicInfo>)(std::vector<BondInfo>)(std::vector<AtomicInfo>)(std::vector<Vector>)
+#define statereferences (clonedatoms)(clonedbonds)(movedatoms)(MovedToVector)
+
+// some defines for all the names, you may use ACTION, STATE and PARAMS
+#define CATEGORY Fill
+#define MENUNAME "fill"
+#define MENUPOSITION 1
+#define ACTIONNAME SuspendInMolecule
+#define TOKEN "suspend-in-molecule"
+
+// finally the information stored in the ActionTrait specialization
+#define DESCRIPTION "\
+suspend the simulation domain with already present molecules in the one \
+selected molecule (must be sphere- or cube-like) such that in the domain the \
+mean density is as specified"
+#define SHORTFORM "u"
Index: src/Actions/FillAction/SuspendInMoleculeAction.hpp
===================================================================
--- src/Actions/FillAction/SuspendInMoleculeAction.hpp	(revision 55736be799388502de9a57ac542d215e5ce0565a)
+++ src/Actions/FillAction/SuspendInMoleculeAction.hpp	(revision 55736be799388502de9a57ac542d215e5ce0565a)
@@ -0,0 +1,22 @@
+/*
+ * SuspendInMoleculeAction.hpp
+ *
+ *  Created on: Sep 05, 2014
+ *      Author: heber
+ */
+
+#ifndef SUSPENDINMOLECULEACTION_HPP
+#define SUSPENDINMOLECULEACTION_HPP
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include "Actions/Action.hpp"
+
+#include "SuspendInMoleculeAction.def"
+#include "Action_impl_header.hpp"
+
+#endif // SUSPENDINMOLECULEACTION_HPP
