Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 1953f90f687a36e6e1f3655d77841bf98a36954c)
+++ src/boundary.cpp	(revision f3278bbc247394ddf549f3044faf52fd11b71543)
@@ -863,4 +863,6 @@
 
             // insert into Filling
+
+            // FIXME: gives completely different results if CopyAtoms[..] used instead of Walker, why???
             *out << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
             Filling->AddAtom(CopyAtoms[Walker->nr]);
@@ -930,4 +932,58 @@
       if (!failflag)
         cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
+      // write temporary envelope
+      if ((DoSingleStepOutput && (mol->TesselStruct->TrianglesOnBoundaryCount % SingleStepWidth == 0))) { // if we have a new triangle and want to output each new triangle configuration
+        TriangleMap::iterator runner = mol->TesselStruct->TrianglesOnBoundary.end();
+        runner--;
+        class BoundaryTriangleSet *triangle = runner->second;
+        if (triangle != NULL) {
+          sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
+          if (DoTecplotOutput) {
+            string NameofTempFile(filename);
+            NameofTempFile.append(NumberName);
+            for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+            NameofTempFile.erase(npos, 1);
+            NameofTempFile.append(TecplotSuffix);
+            *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+            tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+            write_tecplot_file(out, tempstream, mol->TesselStruct, mol, TriangleFilesWritten);
+            tempstream->close();
+            tempstream->flush();
+            delete(tempstream);
+          }
+
+          if (DoRaster3DOutput) {
+            string NameofTempFile(filename);
+            NameofTempFile.append(NumberName);
+            for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+            NameofTempFile.erase(npos, 1);
+            NameofTempFile.append(Raster3DSuffix);
+            *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+            tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+            write_raster3d_file(out, tempstream, mol->TesselStruct, mol);
+    //        // include the current position of the virtual sphere in the temporary raster3d file
+    //        // make the circumsphere's center absolute again
+    //        helper.CopyVector(BaseRay->endpoints[0]->node->node);
+    //        helper.AddVector(BaseRay->endpoints[1]->node->node);
+    //        helper.Scale(0.5);
+    //        (*it)->OptCenter.AddVector(&helper);
+    //        Vector *center = mol->DetermineCenterOfAll(out);
+    //        (*it)->OptCenter.SubtractVector(center);
+    //        delete(center);
+    //        // and add to file plus translucency object
+    //        *tempstream << "# current virtual sphere\n";
+    //        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
+    //        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
+    //          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
+    //          << "\t" << RADIUS << "\t1 0 0\n";
+    //        *tempstream << "9\n  terminating special property\n";
+            tempstream->close();
+            tempstream->flush();
+            delete(tempstream);
+          }
+        }
+        if (DoTecplotOutput || DoRaster3DOutput)
+          TriangleFilesWritten++;
+      }
     } else {
       //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
@@ -941,55 +997,4 @@
       baseline = mol->TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
       flag = false;
-    }
-
-    // write temporary envelope
-    if ((DoSingleStepOutput && (mol->TesselStruct->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
-      class BoundaryTriangleSet *triangle = (mol->TesselStruct->TrianglesOnBoundary.end()--)->second;
-      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
-      if (DoTecplotOutput) {
-        string NameofTempFile(filename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(TecplotSuffix);
-        *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_tecplot_file(out, tempstream, mol->TesselStruct, mol, TriangleFilesWritten);
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
-      }
-
-      if (DoRaster3DOutput) {
-        string NameofTempFile(filename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(Raster3DSuffix);
-        *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_raster3d_file(out, tempstream, mol->TesselStruct, mol);
-//        // include the current position of the virtual sphere in the temporary raster3d file
-//        // make the circumsphere's center absolute again
-//        helper.CopyVector(BaseRay->endpoints[0]->node->node);
-//        helper.AddVector(BaseRay->endpoints[1]->node->node);
-//        helper.Scale(0.5);
-//        (*it)->OptCenter.AddVector(&helper);
-//        Vector *center = mol->DetermineCenterOfAll(out);
-//        (*it)->OptCenter.SubtractVector(center);
-//        delete(center);
-//        // and add to file plus translucency object
-//        *tempstream << "# current virtual sphere\n";
-//        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
-//        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
-//          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
-//          << "\t" << RADIUS << "\t1 0 0\n";
-//        *tempstream << "9\n  terminating special property\n";
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
-      }
-      if (DoTecplotOutput || DoRaster3DOutput)
-        TriangleFilesWritten++;
     }
   }
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 1953f90f687a36e6e1f3655d77841bf98a36954c)
+++ src/boundary.hpp	(revision f3278bbc247394ddf549f3044faf52fd11b71543)
@@ -17,4 +17,5 @@
 #define DEBUG 1
 #define DoSingleStepOutput 0
+#define SingleStepWidth 1
 #define DoTecplotOutput 1
 #define DoRaster3DOutput 1
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 1953f90f687a36e6e1f3655d77841bf98a36954c)
+++ src/builder.cpp	(revision f3278bbc247394ddf549f3044faf52fd11b71543)
@@ -1348,5 +1348,5 @@
             cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
             cout << "\t-b xx xy xz yy yz zz\tCenter atoms in domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
-            cout << "\t-B <basis>\tSetting basis to store to MPQC config files." << endl;
+            cout << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
             cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
             cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
@@ -1358,4 +1358,5 @@
             cout << "\t-h/-H/-?\tGive this help screen." << endl;
             cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
+            cout << "\t-M <basis>\tSetting basis to store to MPQC config files." << endl;
             cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
             cout << "\t-N <radius> <file>\tGet non-convex-envelope." << endl;
@@ -1506,5 +1507,5 @@
         if (config_present == present) {
           switch(argv[argptr-1][1]) {
-            case 'B':
+            case 'M':
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
@@ -1767,4 +1768,21 @@
                 // center
                 mol->CenterInBox((ofstream *)&cout);
+                argptr+=6;
+              }
+              break;
+            case 'B':
+              ExitFlag = 1;
+              if ((argptr+5 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl;
+              } else {
+                SaveFlag = true;
+                j = -1;
+                cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
+                for (int i=0;i<6;i++) {
+                  mol->cell_size[i] = atof(argv[argptr+i]);
+                }
+                // center
+                mol->BoundInBox((ofstream *)&cout);
                 argptr+=6;
               }
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision 1953f90f687a36e6e1f3655d77841bf98a36954c)
+++ src/molecules.cpp	(revision f3278bbc247394ddf549f3044faf52fd11b71543)
@@ -741,4 +741,38 @@
   delete(Minv);
   delete(Center);
+  return status;
+};
+
+
+/** Bounds the molecule in the box whose lengths are defined by vector \a *BoxLengths.
+ * \param *out output stream for debugging
+ */
+bool molecule::BoundInBox(ofstream *out)
+{
+  bool status = true;
+  Vector x;
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *Minv = x.InverseMatrix(M);
+
+  // go through all atoms
+  atom *ptr = start;  // start at first in list
+  while (ptr->next != end) {  // continue with second if present
+    ptr = ptr->next;
+    //ptr->Output(1,1,out);
+    // multiply its vector with matrix inverse
+    x.CopyVector(&ptr->x);
+    x.MatrixMultiplication(Minv);
+    // truncate to [0,1] for each axis
+    for (int i=0;i<NDIM;i++) {
+      while (x.x[i] >= 1.)
+        x.x[i] -= 1.;
+      while (x.x[i] < 0.)
+        x.x[i] += 1.;
+    }
+    x.MatrixMultiplication(M);
+    ptr->x.CopyVector(&x);
+  }
+  delete(M);
+  delete(Minv);
   return status;
 };
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision 1953f90f687a36e6e1f3655d77841bf98a36954c)
+++ src/molecules.hpp	(revision f3278bbc247394ddf549f3044faf52fd11b71543)
@@ -175,4 +175,5 @@
   void CalculateOrbitals(class config &configuration);
   bool CenterInBox(ofstream *out);
+  bool BoundInBox(ofstream *out);
   void CenterEdge(ofstream *out, Vector *max);
   void CenterOrigin(ofstream *out);
