Changeset 36ec71

Timestamp:

Jul 24, 2009, 10:38:32 AM (16 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:

d30402

Parents:

042f82 (diff), 51c910 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

git-author:

Frederik Heber <heber@…> (07/23/09 14:23:32)

git-committer:

Frederik Heber <heber@…> (07/24/09 10:38:32)

Message:

Merge branch 'master' into ConcaveHull

Conflicts:

molecuilder/src/analyzer.cpp
molecuilder/src/bond.cpp
molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/datacreator.cpp
molecuilder/src/datacreator.hpp
molecuilder/src/defs.hpp
molecuilder/src/ellipsoid.cpp
molecuilder/src/joiner.cpp
molecuilder/src/molecules.cpp
molecuilder/src/molecules.hpp
molecuilder/src/parser.cpp
molecuilder/src/parser.hpp

merges:

• analyzer.cpp: all from master, due to Hessian; only Chi(PAS)Fragments from HEAD
• bond.cpp: bond::GetOtherAtom gives cerr output if neither atom matches from master
• boundary.cpp: all from HEAD, is tesselation stuff
• boundary.hpp: all from HEAD, is tesselation stuff
• builder.cpp: Center...(), RemoveAtoms(), SaveConfig() from master, is MultipleMolecules, rest from HEAD (Translation and ...)
• config.cpp: FileBuffer stuff from HEAD, only some taken over from HEAD
• datacreator.cpp: all from master, is Hessian
• datacreator.hpp: all from master, is Hessian
• defs.hpp: all from master, is Thermostat
• ellipsoid.cpp: all from HEAD, created during tesselation
• joiner.cpp: all from master, is Hessian
• molecules.cpp: molecule::molecule(), ConstrainedPotential() stuff, Thermostats(), RemoveAtom(),CreateAdjacencyList(), FragmentMolecule() from master, CenterInBox(), CenterEdge(), CenterOrigin(), CenterPeriodic() (instead of CenterAtVector()), TranslatePeriodically(), DeterminePeriodicCenter(), Output...(), ReturnFullForSymmnetricMatrix() from HEAD
• molecules.hpp: see above, ConfigfileBuffer struct from HEAD, class config joined
• parser.cpp: all from master, is Hessian
• parser.hpp: all from master, is Hessian

compilation fixes:

• molecules.cpp: PrincipalAxisSystem() - CenterGravity() -> CenterOrigin(), DetermineCenter() -> DeterminePeriodicCenter
• datacreator.cpp: all char * to const char *
• builder.cpp: ParseCommandLineOptions() - case 'c': had argptr++ and argptr+=3; case 'O' - x was not zero'd

Location:

src

Files:

• Makefile.am (modified) (1 diff)
• analyzer.cpp (modified) (12 diffs)
• atom.cpp (modified) (1 diff)
• bond.cpp (modified) (2 diffs)
• boundary.cpp (modified) (2 diffs)
• builder.cpp (modified) (37 diffs)
• config.cpp (modified) (10 diffs)
• datacreator.cpp (modified) (37 diffs)
• datacreator.hpp (modified) (1 diff)
• defs.hpp (modified) (1 diff)
• graph.cpp (modified) (1 diff)
• graph.hpp (added)
• joiner.cpp (modified) (10 diffs)
• moleculelist.cpp (modified) (4 diffs)
• molecules.cpp (modified) (20 diffs)
• molecules.hpp (modified) (9 diffs)
• parser.cpp (modified) (23 diffs)
• parser.hpp (modified) (3 diffs)

src/Makefile.am

@@ -10 +10 @@
 
 
-#EXTRA_DIST = ${molecuilder_DATA}
+EXTRA_DIST = ${molecuilder_DATA}

src/analyzer.cpp

@@ -25 +25 @@
   periodentafel *periode = NULL; // and a period table of all elements
   EnergyMatrix Energy;
+  EnergyMatrix EnergyFragments;
+  ForceMatrix Force;
+  ForceMatrix ForceFragments;
+  HessianMatrix Hessian;
+  HessianMatrix HessianFragments;
   EnergyMatrix Hcorrection;
-  ForceMatrix Force;
+  EnergyMatrix HcorrectionFragments;
   ForceMatrix Shielding;
   ForceMatrix ShieldingPAS;
+  ForceMatrix Chi;
+  ForceMatrix ChiPAS;
   EnergyMatrix Time;
-  EnergyMatrix EnergyFragments;
-  EnergyMatrix HcorrectionFragments;
-  ForceMatrix ForceFragments;
   ForceMatrix ShieldingFragments;
   ForceMatrix ShieldingPASFragments;
+  ForceMatrix ChiFragments;
+  ForceMatrix ChiPASFragments;
   KeySetsContainer KeySet;
   ofstream output;
@@ -49 +55 @@
   stringstream yrange;
   char *dir = NULL;
-  bool Hcorrected = true;
+  bool NoHessian = false;
+  bool NoTime = false;
+  bool NoHCorrection = true;
   int counter;
 
@@ -83 +91 @@
   // ------------- Parse through all Fragment subdirs --------
   if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
-  Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
+  if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0)) {
+    NoHCorrection = true;
+    cout << "No HCorrection file found, skipping these." << endl;
+  }
+  
   if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
-  if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
+  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix,0,0)) {
+    NoHessian = true;
+    cout << "No Hessian file found, skipping these." << endl;
+  }
+  if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) {
+    NoTime = true;
+    cout << "No speed file found, skipping these." << endl;
+  }
   if (periode != NULL) { // also look for PAS values
     if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
     if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
+    if (!Chi.ParseFragmentMatrix(argv[1], dir, ChiSuffix, 1, 0)) return 1;
+    if (!ChiPAS.ParseFragmentMatrix(argv[1], dir, ChiPASSuffix, 1, 0)) return 1;
   }
 
   // ---------- Parse the TE Factors into an array -----------------
   if (!Energy.ParseIndices()) return 1;
-  if (Hcorrected) Hcorrection.ParseIndices();
+  if (!NoHCorrection) Hcorrection.ParseIndices();
 
   // ---------- Parse the Force indices into an array ---------------
   if (!Force.ParseIndices(argv[1])) return 1;
   if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
-  if (!ForceFragments.ParseIndices(argv[1])) return 1;
+  if (!ForceFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
+
+  // ---------- Parse hessian indices into an array -----------------
+  if (!NoHessian) {
+    if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
+    if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1;
+    if (!HessianFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
+  }
 
   // ---------- Parse the shielding indices into an array ---------------
@@ -108 +136 @@
     if(!ShieldingFragments.ParseIndices(argv[1])) return 1;
     if(!ShieldingPASFragments.ParseIndices(argv[1])) return 1;
+    if(!Chi.ParseIndices(argv[1])) return 1;
+    if(!ChiPAS.ParseIndices(argv[1])) return 1;
+    if (!ChiFragments.AllocateMatrix(Chi.Header, Chi.MatrixCounter, Chi.RowCounter, Chi.ColumnCounter)) return 1;
+    if (!ChiPASFragments.AllocateMatrix(ChiPAS.Header, ChiPAS.MatrixCounter, ChiPAS.RowCounter, ChiPAS.ColumnCounter)) return 1;
+    if(!ChiFragments.ParseIndices(argv[1])) return 1;
+    if(!ChiPASFragments.ParseIndices(argv[1])) return 1;
   }
 
@@ -116 +150 @@
   // ---------- Parse fragment files created by 'joiner' into an array -------------
   if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
-  if (Hcorrected) HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
+  if (!NoHCorrection) 
+    HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
   if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
+  if (!NoHessian)
+    if (!HessianFragments.ParseFragmentMatrix(argv[1], dir, HessianFragmentSuffix,0,0)) return 1;
   if (periode != NULL) { // also look for PAS values
     if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
     if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASFragmentSuffix, 1, 0)) return 1;
+    if (!ChiFragments.ParseFragmentMatrix(argv[1], dir, ChiFragmentSuffix, 1, 0)) return 1;
+    if (!ChiPASFragments.ParseFragmentMatrix(argv[1], dir, ChiPASFragmentSuffix, 1, 0)) return 1;
   }
 
@@ -129 +168 @@
   filename << argv[3] << "/" << "energy-forces.all";
   output.open(filename.str().c_str(), ios::out);
-  output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header << endl;
+  output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header[Energy.MatrixCounter] << endl;
   for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
-    for(int k=0;k<Energy.ColumnCounter;k++)
+    for(int k=0;k<Energy.ColumnCounter[Energy.MatrixCounter];k++)
       output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
     output << endl;
   }
   output << endl;
-
-  output << endl << "Total Forces" << endl << "===============" << endl << Force.Header << endl;
+  
+  output << endl << "Total Forces" << endl << "===============" << endl << Force.Header[Force.MatrixCounter] << endl;
   for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
-    for(int k=0;k<Force.ColumnCounter;k++)
+    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
       output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
     output << endl;
@@ -145 +184 @@
   output << endl;
 
+  if (!NoHessian) {
+    output << endl << "Total Hessian" << endl << "===============" << endl << Hessian.Header[Hessian.MatrixCounter] << endl;
+    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
+      for(int k=0;k<Hessian.ColumnCounter[Hessian.MatrixCounter];k++)
+        output << scientific << Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output << endl;
+  }
+
   if (periode != NULL) { // also look for PAS values
-    output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header << endl;
+    output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header[Hessian.MatrixCounter] << endl;
     for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
-      for(int k=0;k<Shielding.ColumnCounter;k++)
+      for(int k=0;k<Shielding.ColumnCounter[Shielding.MatrixCounter];k++)
         output << scientific << Shielding.Matrix[ Shielding.MatrixCounter ][j][k] << "\t";
       output << endl;
     }
     output << endl;
-
-    output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header << endl;
+  
+    output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header[ShieldingPAS.MatrixCounter] << endl;
     for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
-      for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
+      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
         output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
       output << endl;
     }
     output << endl;
-  }
-
-  output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
-  for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
-    for(int k=0;k<Time.ColumnCounter;k++) {
-      output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
-    }
-    output << endl;
-  }
-  output << endl;
+
+    output << endl << "Total Chis" << endl << "===============" << endl << Chi.Header[Chi.MatrixCounter] << endl;
+    for(int j=0;j<Chi.RowCounter[Chi.MatrixCounter];j++) {
+      for(int k=0;k<Chi.ColumnCounter[Chi.MatrixCounter];k++)
+        output << scientific << Chi.Matrix[ Chi.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output << endl;
+  
+    output << endl << "Total Chis PAS" << endl << "===============" << endl << ChiPAS.Header[ChiPAS.MatrixCounter] << endl;
+    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
+      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
+        output << scientific << ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output << endl;
+  }
+  
+  if (!NoTime) {
+    output << endl << "Total Times" << endl << "===============" << endl << Time.Header[Time.MatrixCounter] << endl;
+    for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
+      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
+        output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
+      }
+      output << endl;
+    }
+    output << endl;
+  }
   output.close();
-  for(int k=0;k<Time.ColumnCounter;k++)
-    Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
+  if (!NoTime)
+    for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++)
+      Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
 
   // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
@@ -183 +251 @@
   // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
   // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
-  if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
-  if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
-  for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
-    for(int k=Time.ColumnCounter;k--;) {
-      Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
-    }
-  counter = 0;
-  output << "#Order\tFrag.No.\t" << Time.Header << endl;
-  output2 << "#Order\tFrag.No.\t" << Time.Header << endl;
-  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
-    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
-      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
-        for(int k=Time.ColumnCounter;k--;) {
-          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
-        }
-    counter += KeySet.FragmentsPerOrder[BondOrder];
-    output << BondOrder+1 << "\t" << counter;
-    output2 << BondOrder+1 << "\t" << counter;
-    for(int k=0;k<Time.ColumnCounter;k++) {
-      output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
-      if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
-        output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
-      else
-        output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
-    }
-    output << endl;
-    output2 << endl;
-  }
-  output.close();
-  output2.close();
+  if (!NoTime) {
+    if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
+    if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
+    for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
+      for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
+        Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
+      }
+    counter = 0;
+    output << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
+    output2 << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
+    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+      for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
+        for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
+          for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
+            Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+          }
+      counter += KeySet.FragmentsPerOrder[BondOrder];
+      output << BondOrder+1 << "\t" << counter;
+      output2 << BondOrder+1 << "\t" << counter;
+      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
+        output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
+        if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
+          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
+        else
+          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
+      }
+      output << endl;
+      output2 << endl;
+    }
+    output.close();
+    output2.close();
+  }
+
+  if (!NoHessian) {
+    // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in hessian to full QM
+    if (!CreateDataDeltaHessianOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaHessian_xx-Order", "Plot of error between approximated hessian and full hessian versus the Bond Order", datum)) return 1;
+
+    if (!CreateDataDeltaFrobeniusOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaFrobeniusHessian_xx-Order", "Plot of error between approximated hessian and full hessian in the frobenius norm versus the Bond Order", datum)) return 1;
+
+    // ++++++++++++++++++++++++++++++++++++++Plotting Hessian vs. Order
+    if (!CreateDataHessianOrderPerAtom(HessianFragments, KeySet, argv[3], "Hessian_xx-Order", "Plot of approximated hessian versus the Bond Order", datum)) return 1;
+    if (!AppendOutputFile(output, argv[3], "Hessian_xx-Order.dat" )) return false;
+    output << endl << "# Full" << endl;
+    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
+      output << j << "\t";
+      for(int k=0;k<Hessian.ColumnCounter[Force.MatrixCounter];k++)
+        output << scientific <<  Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output.close();
+  }
 
   // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings
-
   if (periode != NULL) { // also look for PAS values
     if (!CreateDataDeltaForcesOrderPerAtom(ShieldingPAS, ShieldingPASFragments, KeySet, argv[3], "DeltaShieldingsPAS-Order", "Plot of error between approximated shieldings and full shieldings versus the Bond Order", datum)) return 1;
@@ -223 +311 @@
     for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
       output << j << "\t";
-      for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
+      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
         output << scientific <<  ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
       output << endl;
     }
-  }
-  output.close();
+    output.close();
+    if (!CreateDataDeltaForcesOrderPerAtom(ChiPAS, ChiPASFragments, KeySet, argv[3], "DeltaChisPAS-Order", "Plot of error between approximated Chis and full Chis versus the Bond Order", datum)) return 1;
+    if (!CreateDataForcesOrderPerAtom(ChiPASFragments, KeySet, argv[3], "ChisPAS-Order", "Plot of approximated Chis versus the Bond Order", datum)) return 1;
+    if (!AppendOutputFile(output, argv[3], "ChisPAS-Order.dat" )) return false;
+    output << endl << "# Full" << endl;
+    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
+      output << j << "\t";
+      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
+        output << scientific <<  ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
+      output << endl;
+    }
+    output.close();
+  }
 
 
@@ -255 +354 @@
   for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
     output << j << "\t";
-    for(int k=0;k<Force.ColumnCounter;k++)
+    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
       output << scientific <<  Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
     output << endl;
@@ -303 +402 @@
   Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
   yrange.str("[1e-8:1e+1]");
-
-  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
-  if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "",  1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
-
+  
+  if (!NoTime) {
+    // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
+    if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "",  1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
+  }
+  
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
   if (!CreatePlotOrder(Energy, KeySet, argv[3], "DeltaEnergies-Order", 1, "outside", "y", "",  1, 1, "bond order k", "absolute error in energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
@@ -403 +504 @@
     }
     output << "'ShieldingsPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
-    output.close();
-    output2.close();
+    output2.close();  
+
+    if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
+    if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
+    CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
+    CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
+    output << "set boxwidth " << step << endl;
+    output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
+    output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
+    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+      output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
+      output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
+      if (BondOrder-1 != KeySet.Order)
+        output2 << ", \\" << endl;
+    }
+    output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
+    output.close();  
+    output2.close();  
+
+    if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
+    if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
+    CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
+    CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
+    output << "set boxwidth " << step << endl;
+    output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
+    output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
+    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+      output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
+      output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
+      if (BondOrder-1 != KeySet.Order)
+        output2 << ", \\" << endl;
+    }
+    output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
+    output.close();  
+    output2.close();  
   }
 

src/atom.cpp

@@ -89 +89 @@
 };
 
-ostream & operator << (ostream &ost, atom &a)
+ostream & operator << (ostream &ost, const atom &a) 
 {
   ost << "[" << a.Name << "|" << &a << "]";

src/bond.cpp

@@ -82 +82 @@
 };
 
-ostream & operator << (ostream &ost, bond &b) 
+ostream & operator << (ostream &ost, const bond &b) 
 {
   ost << "[" << b.leftatom->Name << " <" << b.BondDegree << "(H" << b.HydrogenBond << ")>" << b.rightatom->Name << "]";
@@ -98 +98 @@
   if(rightatom == Atom) 
     return leftatom;
+  cerr << "Bond " << *this << " does not contain atom " << *Atom << "!" << endl;
   return NULL;
 };

src/boundary.cpp

@@ -10 +10 @@
 #define DoSingleStepOutput 0
 #define DoTecplotOutput 1
-#define DoRaster3DOutput 0
+#define DoRaster3DOutput 1
 #define DoVRMLOutput 1
 #define TecplotSuffix ".dat"
@@ -45 +45 @@
 void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
-  cout << Verbose(6) << "Adding line " << *line << " to " << *this << "." << endl;
+  cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
+      << endl;
   if (line->endpoints[0] == this)
     {

src/builder.cpp

@@ -282 +282 @@
     case 'a':
       cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
-      mol->CenterOrigin((ofstream *)&cout, &x);
+      mol->CenterOrigin((ofstream *)&cout);
       break;
     case 'b':
       cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
-      mol->CenterGravity((ofstream *)&cout, &x);
+      mol->CenterPeriodic((ofstream *)&cout);
       break;
     case 'c':
@@ -295 +295 @@
       }
       mol->CenterEdge((ofstream *)&cout, &x);  // make every coordinate positive
-      mol->Translate(&y); // translate by boundary
+      mol->Center.AddVector(&y); // translate by boundary
       helper.CopyVector(&y);
       helper.Scale(2.);
@@ -307 +307 @@
         cin >> x.x[i];
       }
+      // update Box of atoms by boundary
+      mol->SetBoxDimension(&x);
       // center
       mol->CenterInBox((ofstream *)&cout);
-      // update Box of atoms by boundary
-      mol->SetBoxDimension(&x);
       break;
   }
@@ -463 +463 @@
       cin >> tmp1;
       first = mol->start;
-      while(first->next != mol->end) {
-        first = first->next;
+      second = first->next;
+      while(second != mol->end) {
+        first = second;
+        second = first->next;
         if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
           mol->RemoveAtom(first);
@@ -472 +474 @@
       cout << Verbose(0) << "Which axis is it: ";
       cin >> axis;
-      cout << Verbose(0) << "Left inward boundary: ";
+      cout << Verbose(0) << "Lower boundary: ";
       cin >> tmp1;
-      cout << Verbose(0) << "Right inward boundary: ";
+      cout << Verbose(0) << "Upper boundary: ";
       cin >> tmp2;
       first = mol->start;
-      while(first->next != mol->end) {
-        first = first->next;
-        if ((first->x.x[axis] > tmp2) || (first->x.x[axis] < tmp1)) // out of boundary ...
+      second = first->next;
+      while(second != mol->end) {
+        first = second;
+        second = first->next;
+        if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
+          //cout << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
           mol->RemoveAtom(first);
+        }
       }
       break;
@@ -644 +650 @@
   cout << Verbose(0) << "all else - go back" << endl;
   cout << Verbose(0) << "===============================================" << endl;
-  if (molecules->NumberOfActiveMolecules() > 0)
+  if (molecules->NumberOfActiveMolecules() > 1)
     cout << Verbose(0) << "WARNING: There is more than one molecule active! Atoms will be added to each." << endl;
   cout << Verbose(0) << "INPUT: ";
@@ -655 +661 @@
 
     case 'a': // add atom
-      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         mol = *ListRunner;
         cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
@@ -663 +670 @@
 
     case 'b': // scale a bond
-      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         mol = *ListRunner;
         cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
@@ -685 +693 @@
 
     case 'c': // unit scaling of the metric
-      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         mol = *ListRunner;
         cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
@@ -701 +710 @@
 
     case 'l': // measure distances or angles
-      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         mol = *ListRunner;
         cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
@@ -709 +719 @@
 
     case 'r': // remove atom
-      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         mol = *ListRunner;
         cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
@@ -717 +728 @@
 
     case 'u': // change an atom's element
-      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         int Z;
         mol = *ListRunner;
@@ -761 +773 @@
   cout << Verbose(0) << "all else - go back" << endl;
   cout << Verbose(0) << "===============================================" << endl;
-  if (molecules->NumberOfActiveMolecules() > 0)
+  if (molecules->NumberOfActiveMolecules() > 1)
     cout << Verbose(0) << "WARNING: There is more than one molecule active! Atoms will be added to each." << endl;
   cout << Verbose(0) << "INPUT: ";
@@ -772 +784 @@
 
     case 'd': // duplicate the periodic cell along a given axis, given times
-      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) 
+        if ((*ListRunner)->ActiveFlag) {
         mol = *ListRunner;
         cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
@@ -830 +843 @@
 
     case 'g': // center the atoms
-      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         mol = *ListRunner;
         cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
@@ -838 +852 @@
 
     case 'i': // align all atoms
-      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         mol = *ListRunner;
         cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
@@ -846 +861 @@
 
     case 'm': // mirror atoms along a given axis
-      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         mol = *ListRunner;
         cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
@@ -854 +870 @@
 
     case 'o': // create the connection matrix
-      {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         double bonddistance;
         clock_t start,end;
@@ -868 +885 @@
 
     case 't': // translate all atoms
-     for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if ((*ListRunner)->ActiveFlag) {
         mol = *ListRunner;
         cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
         cout << Verbose(0) << "Enter translation vector." << endl;
         x.AskPosition(mol->cell_size,0);
-        mol->Translate((const Vector *)&x);
+        mol->Center.AddVector((const Vector *)&x);
      }
      break;
@@ -895 +913 @@
 {
   char choice;  // menu choice char
-  Vector Center;
+  Vector center;
   int nr, count;
   molecule *mol = NULL;
-  char filename[MAXSTRINGSIZE];
-
-  cout << Verbose(0) << "==========Edit MOLECULES=====================" << endl;
+
+  cout << Verbose(0) << "==========EDIT MOLECULES=====================" << endl;
   cout << Verbose(0) << "c - create new molecule" << endl;
   cout << Verbose(0) << "l - load molecule from xyz file" << endl;
   cout << Verbose(0) << "n - change molecule's name" << endl;
   cout << Verbose(0) << "N - give molecules filename" << endl;
-  cout << Verbose(0) << "p - parse xyz file into molecule" << endl;
+  cout << Verbose(0) << "p - parse atoms in xyz file into molecule" << endl;
   cout << Verbose(0) << "r - remove a molecule" << endl;
   cout << Verbose(0) << "all else - go back" << endl;
@@ -921 +938 @@
       break;
 
-    case 'l': // laod from XYZ file
-      cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
-      mol = new molecule(periode);
-      do {
-        cout << Verbose(0) << "Enter file name: ";
+    case 'l': // load from XYZ file
+      {
+        char filename[MAXSTRINGSIZE];
+        cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
+        mol = new molecule(periode);
+        do {
+          cout << Verbose(0) << "Enter file name: ";
+          cin >> filename;
+        } while (!mol->AddXYZFile(filename));
+        mol->SetNameFromFilename(filename);
+        // center at set box dimensions
+        mol->CenterEdge((ofstream *)&cout, &center);
+        mol->cell_size[0] = center.x[0];
+        mol->cell_size[1] = 0;
+        mol->cell_size[2] = center.x[1];
+        mol->cell_size[3] = 0;
+        mol->cell_size[4] = 0;
+        mol->cell_size[5] = center.x[2];
+        molecules->insert(mol);
+      }
+      break;
+
+    case 'n':
+      {
+        char filename[MAXSTRINGSIZE];
+        do {
+          cout << Verbose(0) << "Enter index of molecule: ";
+          cin >> nr;
+          mol = molecules->ReturnIndex(nr);
+        } while (mol == NULL);
+        cout << Verbose(0) << "Enter name: ";
         cin >> filename;
-      } while (!mol->AddXYZFile(filename));
-      mol->SetNameFromFilename(filename);
-      // center at set box dimensions
-      mol->CenterEdge((ofstream *)&cout, &Center);
-      mol->cell_size[0] = Center.x[0];
-      mol->cell_size[1] = 0;
-      mol->cell_size[2] = Center.x[1];
-      mol->cell_size[3] = 0;
-      mol->cell_size[4] = 0;
-      mol->cell_size[5] = Center.x[2];
-      molecules->insert(mol);
-      break;
-
-    case 'n':
-      do {
-        cout << Verbose(0) << "Enter index of molecule: ";
-        cin >> nr;
-        mol = molecules->ReturnIndex(nr);
-      } while (mol != NULL);
-      cout << Verbose(0) << "Enter name: ";
-      cin >> filename;
-      strcpy(mol->name, filename);
+        strcpy(mol->name, filename);
+      }
       break;
 
     case 'N':
-      do {
-        cout << Verbose(0) << "Enter index of molecule: ";
-        cin >> nr;
-        mol = molecules->ReturnIndex(nr);
-      } while (mol != NULL);
-      cout << Verbose(0) << "Enter name: ";
-      cin >> filename;
-      mol->SetNameFromFilename(filename);
+      {
+        char filename[MAXSTRINGSIZE];
+        do {
+          cout << Verbose(0) << "Enter index of molecule: ";
+          cin >> nr;
+          mol = molecules->ReturnIndex(nr);
+        } while (mol == NULL);
+        cout << Verbose(0) << "Enter name: ";
+        cin >> filename;
+        mol->SetNameFromFilename(filename);
+      }
       break;
 
     case 'p': // parse XYZ file
-      mol = NULL;
-      do {
-        cout << Verbose(0) << "Enter index of molecule: ";
-        cin >> nr;
-        mol = molecules->ReturnIndex(nr);
-      } while (mol == NULL);
-      cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
-      do {
-        cout << Verbose(0) << "Enter file name: ";
-        cin >> filename;
-      } while (!mol->AddXYZFile(filename));
-      mol->SetNameFromFilename(filename);
+      {
+        char filename[MAXSTRINGSIZE];
+        mol = NULL;
+        do {
+          cout << Verbose(0) << "Enter index of molecule: ";
+          cin >> nr;
+          mol = molecules->ReturnIndex(nr);
+        } while (mol == NULL);
+        cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
+        do {
+          cout << Verbose(0) << "Enter file name: ";
+          cin >> filename;
+        } while (!mol->AddXYZFile(filename));
+        mol->SetNameFromFilename(filename);
+      }
       break;
 
@@ -981 +1010 @@
       cin >> nr;
       count = 1;
-      MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin();
-      for(; ((ListRunner != molecules->ListOfMolecules.end()) && (count < nr)); ListRunner++);
-      mol = *ListRunner;
-      if (count == nr) {
-        molecules->ListOfMolecules.erase(ListRunner);
-        delete(mol);
-      }
+      for( MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+        if (nr == (*ListRunner)->IndexNr) {
+          mol = *ListRunner;
+          molecules->ListOfMolecules.erase(ListRunner);
+          delete(mol);
+        }
       break;
   }
@@ -1002 +1030 @@
 
   cout << Verbose(0) << "===========MERGE MOLECULES=====================" << endl;
+  cout << Verbose(0) << "a - simple add of one molecule to another" << endl;
   cout << Verbose(0) << "e - embedding merge of two molecules" << endl;
   cout << Verbose(0) << "m - multi-merge of all molecules" << endl;
@@ -1016 +1045 @@
       break;
 
+    case 'a':
+      {
+        int src, dest;
+        molecule *srcmol = NULL, *destmol = NULL;
+        {
+          do {
+            cout << Verbose(0) << "Enter index of destination molecule: ";
+            cin >> dest;
+            destmol = molecules->ReturnIndex(dest);
+          } while ((destmol == NULL) && (dest != -1));
+          do {
+            cout << Verbose(0) << "Enter index of source molecule to add from: ";
+            cin >> src;
+            srcmol = molecules->ReturnIndex(src);
+          } while ((srcmol == NULL) && (src != -1));
+          if ((src != -1) && (dest != -1))
+            molecules->SimpleAdd(srcmol, destmol);
+        }
+      }
+      break;
+
     case 'e':
+      cout << Verbose(0) << "Not implemented yet." << endl;
       break;
 
     case 'm':
+      {
+        int nr;
+        molecule *mol = NULL;
+        do {
+          cout << Verbose(0) << "Enter index of molecule to merge into: ";
+          cin >> nr;
+          mol = molecules->ReturnIndex(nr);
+        } while ((mol == NULL) && (nr != -1));
+        if (nr != -1) {
+          int N = molecules->ListOfMolecules.size()-1;
+          int *src = new int(N);
+          for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+            if ((*ListRunner)->IndexNr != nr)
+              src[N++] = (*ListRunner)->IndexNr;        
+          molecules->SimpleMultiMerge(mol, src, N);
+          delete[](src);
+        }
+      }
       break;
 
     case 's':
+      cout << Verbose(0) << "Not implemented yet." << endl;
       break;
 
     case 't':
+      {
+        int src, dest;
+        molecule *srcmol = NULL, *destmol = NULL;
+        {
+          do {
+            cout << Verbose(0) << "Enter index of destination molecule: ";
+            cin >> dest;
+            destmol = molecules->ReturnIndex(dest);
+          } while ((destmol == NULL) && (dest != -1));
+          do {
+            cout << Verbose(0) << "Enter index of source molecule to merge into: ";
+            cin >> src;
+            srcmol = molecules->ReturnIndex(src);
+          } while ((srcmol == NULL) && (src != -1));
+          if ((src != -1) && (dest != -1))
+            molecules->SimpleMerge(srcmol, destmol);
+        }
+      }
       break;
   }
@@ -1125 +1213 @@
   molecule *mol = new molecule(periode);
 
-  // merge all molecules in MoleculeListClass into this molecule
+  // translate each to its center and merge all molecules in MoleculeListClass into this molecule
   int N = molecules->ListOfMolecules.size();
   int *src = new int(N);
   N=0;
-  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     src[N++] = (*ListRunner)->IndexNr;
+    (*ListRunner)->Translate(&(*ListRunner)->Center);
+  }
   molecules->SimpleMultiAdd(mol, src, N);
+  delete[](src);
+  // ... and translate back
+  for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+    (*ListRunner)->Center.Scale(-1.);
+    (*ListRunner)->Translate(&(*ListRunner)->Center);
+    (*ListRunner)->Center.Scale(-1.);
+  }
 
   cout << Verbose(0) << "Storing configuration ... " << endl;
@@ -1137 +1234 @@
   mol->CalculateOrbitals(*configuration);
   configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
-  strcpy(filename, ConfigFileName);
   if (ConfigFileName != NULL) { // test the file name
-    output.open(ConfigFileName, ios::trunc);
+    strcpy(filename, ConfigFileName);
+    output.open(filename, ios::trunc);
   } else if (strlen(configuration->configname) != 0) {
     strcpy(filename, configuration->configname);
@@ -1227 +1324 @@
   // simply create a new molecule, wherein the config file is loaded and the manipulation takes place
   molecule *mol = new molecule(periode);
+  mol->ActiveFlag = true;
   molecules->insert(mol);
 
@@ -1257 +1355 @@
             cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
             cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
-            cout << "\t-N\tGet non-convex-envelope." << endl;
+            cout << "\t-N <radius> <file>\tGet non-convex-envelope." << endl;
             cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
             cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
+            cout << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl; 
             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
+            cout << "\t-R\t\tRemove all atoms out of sphere around a given one." << endl;
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
             cout << "\t-T x1 x2 x3\tTranslate periodically all atoms by this vector (x1,x2,x3)." << endl;
@@ -1504 +1604 @@
               }
               break;
+            case 'L':
+              ExitFlag = 1;
+              SaveFlag = true;
+              cout << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;
+              if (!mol->LinearInterpolationBetweenConfiguration((ofstream *)&cout, atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration))
+                cout << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;
+              else
+                cout << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;
+              argptr+=3;
+              break;
             case 'P':
               ExitFlag = 1;
@@ -1512 +1622 @@
                 SaveFlag = true;
                 cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
-                if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
+                if (!mol->VerletForceIntegration((ofstream *)&cout, argv[argptr], configuration))
                   cout << Verbose(2) << "File not found." << endl;
                 else
                   cout << Verbose(2) << "File found and parsed." << endl;
                 argptr+=1;
+              }
+              break;
+            case 'R':
+              ExitFlag = 1;
+              if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])))  {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for removing atoms: -R <id> <distance>" << endl;
+              } else {
+                SaveFlag = true;
+                cout << Verbose(1) << "Removing atoms around " << argv[argptr] << " with radius " << argv[argptr+1] << "." << endl;
+                double tmp1 = atof(argv[argptr+1]);
+                atom *third = mol->FindAtom(atoi(argv[argptr]));
+                atom *first = mol->start;
+                if ((third != NULL) && (first != mol->end)) {
+                  atom *second = first->next;
+                  while(second != mol->end) {
+                    first = second;
+                    second = first->next;
+                    if (first->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ...
+                      mol->RemoveAtom(first);
+                  }
+                } else {
+                  cerr << "Removal failed due to missing atoms on molecule or wrong id." << endl;
+                }
+                argptr+=2;
               }
               break;
@@ -1533 +1668 @@
                 argptr+=3;
               }
-              break;
             case 'T':
               ExitFlag = 1;
@@ -1583 +1717 @@
               } 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++) {
@@ -1609 +1744 @@
                 for (int i=0;i<NDIM;i++) {
                   j += i+1;
-                  x.x[i] = atof(argv[argptr++]);
+                  x.x[i] = atof(argv[argptr+i]);
                   mol->cell_size[j] += x.x[i]*2.;
                 }
@@ -1619 +1754 @@
               ExitFlag = 1;
               SaveFlag = true;
-              cout << Verbose(1) << "Centering atoms in origin." << endl;
-              mol->CenterOrigin((ofstream *)&cout, &x);
+              cout << Verbose(1) << "Centering atoms on edge and setting box dimensions." << endl;
+              x.Zero();
+              mol->CenterEdge((ofstream *)&cout, &x);
               mol->SetBoxDimension(&x);
               argptr+=0;
@@ -1802 +1938 @@
   string line;
   char *ConfigFileName = NULL;
-  int j, count;
+  int j;
 
   // =========================== PARSE COMMAND LINE OPTIONS ====================================
@@ -1844 +1980 @@
     cout << Verbose(0) << "============Menu===============================" << endl;
     cout << Verbose(0) << "a - set molecule (in)active" << endl;
-    cout << Verbose(0) << "e - edit new molecules" << endl;
+    cout << Verbose(0) << "e - edit molecules (load, parse, save)" << endl;
     cout << Verbose(0) << "g - globally manipulate atoms in molecule" << endl;
     cout << Verbose(0) << "M - Merge molecules" << endl;
@@ -1863 +1999 @@
           cout << "Enter index of molecule: ";
           cin >> j;
-          count = 1;
-          MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin();
-          for(; ((ListRunner != molecules->ListOfMolecules.end()) && (count < j)); ListRunner++);
-          if (count == j)
-            (*ListRunner)->ActiveFlag = !(*ListRunner)->ActiveFlag;
+          for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+            if ((*ListRunner)->IndexNr == j)
+              (*ListRunner)->ActiveFlag = !(*ListRunner)->ActiveFlag;
         }
         break;

src/config.cpp

@@ -169 +169 @@
   configpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: configpath");
   configname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: configname");
+  ThermostatImplemented = (int *) Malloc((MaxThermostats)*(sizeof(int)), "config constructor: *ThermostatImplemented");
+  ThermostatNames = (char **) Malloc((MaxThermostats)*(sizeof(char *)), "config constructor: *ThermostatNames");
+  for (int j=0;j<MaxThermostats;j++)
+    ThermostatNames[j] = (char *) MallocString(12*(sizeof(char)), "config constructor: ThermostatNames[]");
+  Thermostat = 4;
+  alpha = 0.;
+  ScaleTempStep = 25;
+  TempFrequency = 2.5;
   strcpy(mainname,"pcp");
   strcpy(defaultpath,"not specified");
@@ -175 +183 @@
   configname[0]='\0';
   basis = "3-21G";
+
+  strcpy(ThermostatNames[0],"None");
+  ThermostatImplemented[0] = 1;
+  strcpy(ThermostatNames[1],"Woodcock"); 
+  ThermostatImplemented[1] = 1;
+  strcpy(ThermostatNames[2],"Gaussian"); 
+  ThermostatImplemented[2] = 1;
+  strcpy(ThermostatNames[3],"Langevin"); 
+  ThermostatImplemented[3] = 1;
+  strcpy(ThermostatNames[4],"Berendsen"); 
+  ThermostatImplemented[4] = 1;
+  strcpy(ThermostatNames[5],"NoseHoover"); 
+  ThermostatImplemented[5] = 1;
 
   FastParsing = false;
@@ -187 +208 @@
   DoFullCurrent=0;
   DoWannier=0;
+  DoConstrainedMD=0;
   CommonWannier=0;
   SawtoothStart=0.01;
@@ -195 +217 @@
 
   MaxOuterStep=0;
-  Deltat=1;
+  Deltat=0.01;
   OutVisStep=10;
   OutSrcStep=5;
@@ -247 +269 @@
   Free((void **)&configpath, "config::~config: *configpath");
   Free((void **)&configname, "config::~config: *configname");
+  Free((void **)&ThermostatImplemented, "config::~config: *ThermostatImplemented");
+  for (int j=0;j<MaxThermostats;j++)
+    Free((void **)&ThermostatNames[j], "config::~config: *ThermostatNames[]");
+  Free((void **)&ThermostatNames, "config::~config: **ThermostatNames");
 };
+
+/** Readin of Thermostat related values from parameter file.
+ * \param *source parameter file
+ */
+void config::InitThermostats(ifstream *source)
+{
+  char *thermo = MallocString(12, "IonsInitRead: thermo");
+  int verbose = 0;
+
+  // read desired Thermostat from file along with needed additional parameters
+  if (ParseForParameter(verbose,source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
+    if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
+      if (ThermostatImplemented[0] == 1) {
+        Thermostat = None;
+      } else {
+        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
+        Thermostat = None;
+      }
+    } else if (strcmp(thermo, ThermostatNames[1]) == 0) { // Woodcock
+      if (ThermostatImplemented[1] == 1) {
+        Thermostat = Woodcock;
+        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
+      } else {
+        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
+        Thermostat = None;
+      }
+    } else if (strcmp(thermo, ThermostatNames[2]) == 0) { // Gaussian
+      if (ThermostatImplemented[2] == 1) {
+        Thermostat = Gaussian;
+        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
+      } else {
+        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
+        Thermostat = None;
+      }
+    } else if (strcmp(thermo, ThermostatNames[3]) == 0) { // Langevin
+      if (ThermostatImplemented[3] == 1) {
+        Thermostat = Langevin;
+        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
+        if (ParseForParameter(verbose,source,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
+          cout << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl;
+        } else {
+          alpha = 1.;
+        }
+      } else {
+        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
+        Thermostat = None;
+      }
+    } else if (strcmp(thermo, ThermostatNames[4]) == 0) { // Berendsen
+      if (ThermostatImplemented[4] == 1) {
+        Thermostat = Berendsen;
+        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
+      } else {
+        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
+        Thermostat = None;
+      }
+    } else if (strcmp(thermo, ThermostatNames[5]) == 0) { // Nose-Hoover
+      if (ThermostatImplemented[5] == 1) {
+        Thermostat = NoseHoover;
+        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
+        alpha = 0.;
+      } else {
+        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
+        Thermostat = None;
+      }
+    } else {
+      cout << Verbose(1) << " Warning: thermostat name was not understood!" << endl;
+      Thermostat = None;
+    }
+  } else {
+    if ((MaxOuterStep > 0) && (TargetTemp != 0)) 
+      cout << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl;
+    Thermostat = None;
+  }
+  Free((void **)&thermo, "InitThermostats: thermo");
+};
+
 
 /** Displays menu for editing each entry of the config file.
@@ -597 +699 @@
 void config::Load(char *filename, periodentafel *periode, molecule *mol)
 {
+  ifstream *file = new ifstream(filename);
+  if (file == NULL) {
+    cerr << "ERROR: config file " << filename << " missing!" << endl;
+    return;
+  }
   RetrieveConfigPathAndName(filename);
 
@@ -614 +721 @@
   int verbose = 0;
   double value[3];
-
+  
+  InitThermostats(file);
+  
   /* Namen einlesen */
 
@@ -667 +776 @@
     if (config::SawtoothStart > 1.) config::SawtoothStart = 1.;
   }
-
+  
+  if (ParseForParameter(verbose,FileBuffer,"DoConstrainedMD", 0, 1, 1, int_type, &(config::DoConstrainedMD), 1, optional))
+    if (config::DoConstrainedMD < 0) 
+      config::DoConstrainedMD = 0;
   ParseForParameter(verbose,FileBuffer,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
   if (!ParseForParameter(verbose,FileBuffer,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
@@ -1168 +1280 @@
     *output << "DoPerturbation\t" << config::DoPerturbation << "\t# Do perturbation calculate and determine susceptibility and shielding" << endl;
     *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
+    *output << "DoConstrainedMD\t" << config::DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
+    *output << "Thermostat\t" << ThermostatNames[Thermostat] << "\t";
+    switch(Thermostat) {
+      default:
+      case None:
+        break;
+      case Woodcock:
+        *output << ScaleTempStep;
+        break;
+      case Gaussian:
+        *output << ScaleTempStep;
+        break;
+      case Langevin:
+        *output << TempFrequency << "\t" << alpha;
+        break;
+      case Berendsen:
+        *output << TempFrequency;
+        break;
+      case NoseHoover:
+        *output << HooverMass;
+        break;
+    };
+    *output << "\t# Which Thermostat and its parameters to use in MD case." << endl;
     *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
     *output << "SawtoothStart\t" << config::SawtoothStart << "\t# Absolute value for smooth transition at cell border " << endl;
@@ -1305 +1440 @@
   *output << ")" << endl;
   *output << "basis<GaussianBasisSet>: (" << endl;
-  *output << "\tname = \""<< basis << "\"" << endl;
+  *output << "\tname = \"" << basis << "\"" << endl;
   *output << "\tmolecule = $:molecule" << endl;
   *output << ")" << endl;

src/datacreator.cpp

@@ -63 +63 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
+  output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
       for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
-        for(int k=Fragments.ColumnCounter;k--;)
+        for(int k=Fragments.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];k--;)
           Fragments.Matrix[Fragments.MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     }
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
-    for (int l=0;l<Fragments.ColumnCounter;l++)
+    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
     output << endl;
@@ -96 +96 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
+  output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
       for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
-        for(int k=Fragments.ColumnCounter;k--;)
+        for(int k=Fragments.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];k--;)
           Fragments.Matrix[Fragments.MatrixCounter][j][k] -= Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     }
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
-    for (int l=0;l<Fragments.ColumnCounter;l++)
+    for (int l=0;l<Fragments.ColumnCounter[Energy.MatrixCounter];l++)
       if (fabs(Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]) < MYEPSILON)
         output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
@@ -133 +133 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
+  output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   Fragments.SetLastMatrix(0.,0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
@@ -139 +139 @@
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     CreateForce(Fragments, Fragments.MatrixCounter);
-    for (int l=0;l<Fragments.ColumnCounter;l++)
+    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
        output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
     output << endl;
@@ -165 +165 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
+  output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
@@ -171 +171 @@
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     CreateForce(Fragments, Fragments.MatrixCounter);
-    for (int l=0;l<Fragments.ColumnCounter;l++)
+    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
        output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
     output << endl;
@@ -198 +198 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "# AtomNo\t" << Fragments.Header << endl;
+  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
@@ -207 +207 @@
     for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
       output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
-      for (int l=0;l<Fragments.ColumnCounter;l++) {
+      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
         if (((l+1) % 3) == 0) {
           norm = 0.;
@@ -213 +213 @@
             norm += Force.Matrix[Force.MatrixCounter][ j ][l+m]*Force.Matrix[Force.MatrixCounter][ j ][l+m];
           norm = sqrt(norm);
-        }
+        }                                                                                                    
 //        if (norm < MYEPSILON)
           output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
@@ -244 +244 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "# AtomNo\t" << Fragments.Header << endl;
+  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     //cout << "Current order is " << BondOrder << "." << endl;
@@ -252 +252 @@
     for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
       output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
-      for (int l=0;l<Fragments.ColumnCounter;l++)
+      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++)
         output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
       output << endl;
@@ -262 +262 @@
 };
 
+
+/** Plot hessian error vs. vs atom vs. order.
+ * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
+ * \param &Fragments HessianMatrix class containing matrix values
+ * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param *prefix prefix in filename (without ending)
+ * \param *msg message to be place in first line as a comment
+ * \param *datum current date and time
+ * \return true if file was written successfully
+ */
+bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum) 
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
+  cout << msg << endl;
+  output << "# " << msg << ", created on " << datum;
+  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
+  Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    //cout << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubHessians(Fragments, KeySet, BondOrder, -1.);
+    // errors per atom
+    output << endl << "#Order\t" << BondOrder+1 << endl;
+    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
+      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
+      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
+        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
+      }
+      output << endl;
+    }
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot hessian error vs. vs atom vs. order in the frobenius norm.
+ * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
+ * \param &Fragments HessianMatrix class containing matrix values
+ * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param *prefix prefix in filename (without ending)
+ * \param *msg message to be place in first line as a comment
+ * \param *datum current date and time
+ * \return true if file was written successfully
+ */
+bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum) 
+{
+  stringstream filename;
+  ofstream output;
+  double norm = 0;
+  double tmp;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
+  cout << msg << endl;
+  output << "# " << msg << ", created on " << datum;
+  output << "# AtomNo\t";
+  Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    output << "Order" << BondOrder+1 << "\t";
+  }
+  output << endl;
+  output << Fragments.RowCounter[ Fragments.MatrixCounter ] << "\t";
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    //cout << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubHessians(Fragments, KeySet, BondOrder, -1.);
+    // frobenius norm of errors per atom
+    norm = 0.;
+    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
+      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
+        tmp = Fragments.Matrix[Fragments.MatrixCounter][ j ][l];
+        norm += tmp*tmp;
+      }
+    }
+    output << scientific << sqrt(norm)/(Fragments.RowCounter[ Fragments.MatrixCounter ]*Fragments.ColumnCounter[ Fragments.MatrixCounter] ) << "\t";
+  }
+  output << endl;
+  output.close();
+  return true;
+};
+
+/** Plot hessian error vs. vs atom vs. order.
+ * \param &Fragments HessianMatrix class containing matrix values
+ * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param *prefix prefix in filename (without ending)
+ * \param *msg message to be place in first line as a comment
+ * \param *datum current date and time
+ * \return true if file was written successfully
+ */
+bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum) 
+{
+  stringstream filename;
+  ofstream output;
+
+  filename << prefix << ".dat";
+  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 
+  cout << msg << endl;
+  output << "# " << msg << ", created on " << datum;
+  output << "# AtomNo\t" << Fragments.Header[ Fragments.MatrixCounter ] << endl;
+  Fragments.SetLastMatrix(0., 0);
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    //cout << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubHessians(Fragments, KeySet, BondOrder, 1.);
+    // errors per atom
+    output << endl << "#Order\t" << BondOrder+1 << endl;
+    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 
+      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
+      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++)
+        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
+      output << endl;
+    }
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
 /** Plot matrix vs. fragment.
  */
@@ -273 +393 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum << endl;
-  output << "#Order\tFrag.No.\t" << Fragment.Header << endl;
+  output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
       output << BondOrder+1 << "\t" << KeySet.OrderSet[BondOrder][i]+1;
       CreateFragment(Fragment, KeySet.OrderSet[BondOrder][i]);
-      for (int l=0;l<Fragment.ColumnCounter;l++)
+      for (int l=0;l<Fragment.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];l++)
         output << scientific << "\t" << Fragment.Matrix[ KeySet.OrderSet[BondOrder][i] ][ Fragment.RowCounter[ KeySet.OrderSet[BondOrder][i] ] ][l];
       output << endl;
@@ -297 +417 @@
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
       if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
-        for (int k=Fragments.ColumnCounter;k--;)
+        for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
           Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
       }
@@ -314 +434 @@
     int i=0;
     do {  // first get a minimum value unequal to 0
-      for (int k=Fragments.ColumnCounter;k--;)
+      for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
         Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
       i++;
@@ -320 +440 @@
     for(;i<KeySet.FragmentsPerOrder[BondOrder];i++) { // then find lowest
       if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) > fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
-        for (int k=Fragments.ColumnCounter;k--;)
+        for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
           Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
       }
@@ -338 +458 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "#Order\tFrag.No.\t" << Fragment.Header << endl;
+  output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
   // max
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
@@ -344 +464 @@
     CreateFragmentOrder(Fragment, KeySet, BondOrder);
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
-    for (int l=0;l<Fragment.ColumnCounter;l++)
+    for (int l=0;l<Fragment.ColumnCounter[ Fragment.MatrixCounter ];l++)
       output << scientific << "\t" << Fragment.Matrix[ Fragment.MatrixCounter ][ Fragment.RowCounter[ Fragment.MatrixCounter ]-1 ][l];
     output << endl;
@@ -358 +478 @@
 void CreateEnergy(class MatrixContainer &Energy, int MatrixNumber)
 {
-  for(int k=0;k<Energy.ColumnCounter;k++)
+  for(int k=0;k<Energy.ColumnCounter[MatrixNumber];k++)
     Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber] ] [k] =  Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber]-1 ] [k];
 };
@@ -369 +489 @@
 void CreateMinimumForce(class MatrixContainer &Force, int MatrixNumber)
 {
-  for (int l=Force.ColumnCounter;l--;)
+  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
     Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
-  for (int l=5;l<Force.ColumnCounter;l+=3) {
+  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
     double stored = 0;
     int k=0;
@@ -404 +524 @@
 {
   int divisor = 0;
-  for (int l=Force.ColumnCounter;l--;)
+  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
     Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
-  for (int l=5;l<Force.ColumnCounter;l+=3) {
+  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
     double tmp = 0;
     for (int k=Force.RowCounter[MatrixNumber];k--;) {
@@ -428 +548 @@
 void CreateMaximumForce(class MatrixContainer &Force, int MatrixNumber)
 {
-  for (int l=5;l<Force.ColumnCounter;l+=3) {
+  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
     double stored = 0;
     for (int k=Force.RowCounter[MatrixNumber];k--;) {
@@ -460 +580 @@
 void CreateVectorSumForce(class MatrixContainer &Force, int MatrixNumber)
 {
-  for (int l=Force.ColumnCounter;l--;)
+  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
     Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
-  for (int l=0;l<Force.ColumnCounter;l++) {
+  for (int l=0;l<Force.ColumnCounter[MatrixNumber];l++) {
     for (int k=Force.RowCounter[MatrixNumber];k--;)
       Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] += Force.Matrix[MatrixNumber][k][l];
@@ -538 +658 @@
 void AbsEnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
 {
-  stringstream line(Energy.Header);
+  stringstream line(Energy.Header[ Energy.MatrixCounter ]);
   string token;
 
   getline(line, token, '\t');
-  for (int i=2; i<= Energy.ColumnCounter;i++) {
+  for (int i=2; i<= Energy.ColumnCounter[Energy.MatrixCounter];i++) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
       token.erase(0,1);
     output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(abs($" << i+2 << ")) " << uses;
-    if (i != (Energy.ColumnCounter))
+    if (i != (Energy.ColumnCounter[Energy.MatrixCounter]))
       output << ", \\";
     output << endl;
@@ -562 +682 @@
 void EnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
 {
-  stringstream line(Energy.Header);
+  stringstream line(Energy.Header[Energy.MatrixCounter]);
   string token;
 
   getline(line, token, '\t');
-  for (int i=1; i<= Energy.ColumnCounter;i++) {
+  for (int i=1; i<= Energy.ColumnCounter[Energy.MatrixCounter];i++) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
       token.erase(0,1);
     output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":" << i+2 << " " << uses;
-    if (i != (Energy.ColumnCounter))
+    if (i != (Energy.ColumnCounter[Energy.MatrixCounter]))
       output << ", \\";
     output << endl;
@@ -586 +706 @@
 void ForceMagnitudePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
 {
-  stringstream line(Force.Header);
+  stringstream line(Force.Header[Force.MatrixCounter]);
   string token;
 
@@ -594 +714 @@
   getline(line, token, '\t');
   getline(line, token, '\t');
-  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
@@ -600 +720 @@
     token.erase(token.length(), 1);  // kill residual index char (the '0')
     output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses;
-    if (i != (Force.ColumnCounter-1))
+    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
       output << ", \\";
     output << endl;
@@ -617 +737 @@
 void AbsFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
 {
-  stringstream line(Force.Header);
+  stringstream line(Force.Header[Force.MatrixCounter]);
   string token;
 
@@ -625 +745 @@
   getline(line, token, '\t');
   getline(line, token, '\t');
-  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
@@ -631 +751 @@
     token.erase(token.length(), 1);  // kill residual index char (the '0')
     output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
-    if (i != (Force.ColumnCounter-1))
+    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
       output << ", \\";
     output << endl;
@@ -648 +768 @@
 void BoxesForcePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
 {
-  stringstream line(Force.Header);
-  const char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
+  stringstream line(Force.Header[Force.MatrixCounter]);
+  char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
   string token;
 
@@ -657 +777 @@
   getline(line, token, '\t');
   getline(line, token, '\t');
-  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
@@ -663 +783 @@
     token.erase(token.length(), 1);  // kill residual index char (the '0')
     output << "'" << prefix << ".dat' title '" << token << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses << " " << fillcolor[(i-7)/3];
-    if (i != (Force.ColumnCounter-1))
+    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
       output << ", \\";
     output << endl;
@@ -680 +800 @@
 void BoxesFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
 {
-  stringstream line(Force.Header);
-  const char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
+  stringstream line(Force.Header[Force.MatrixCounter]);
+  char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
   string token;
 
@@ -689 +809 @@
   getline(line, token, '\t');
   getline(line, token, '\t');
-  for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
+  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     getline(line, token, '\t');
     while (token[0] == ' ') // remove leading white spaces
@@ -695 +815 @@
     token.erase(token.length(), 1);  // kill residual index char (the '0')
     output << "'" << prefix << ".dat' title '" << token << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):" << i+1 << " " << uses << " " << fillcolor[(i-7)/3];
-    if (i != (Force.ColumnCounter-1))
+    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
       output << ", \\";
     output << endl;

src/datacreator.hpp

@@ -25 +25 @@
 bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum);
 bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int));
-bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum);
+bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 
+bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 
+bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 
+bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 
 bool CreateDataFragment(class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int));
 bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int));

src/defs.hpp

@@ -16 +16 @@
 #define BONDTHRESHOLD 0.5   //!< CSD threshold in bond check which is the width of the interval whose center is the sum of the covalent radii 
 #define AtomicEnergyToKelvin 315774.67  //!< conversion factor from atomic energy to kelvin via boltzmann factor
+#define KelvinToAtomicTemperature 3.1668152e-06    //!< conversion factor for Kelvin to atomic temperature (Hartree over k_B)
+#define KelvinToeV 8.6173422e-05                   //!< conversion factor for Kelvin to Ht (k_B * T = energy), thus Boltzmann constant in eV/K
+#define AtomicMassUnitsToeV 931494088.        //!< conversion factor for atomic weight in units to mass in eV
+#define AtomicMassUnitsToHt 34480864.        //!< conversion factor for atomic weight in units to mass in Ht (protonmass/electronmass * electron_mass_in_Ht
+#define ElectronMass_Ht 18778.865            //!< electron mass in Ht
+#define ElectronMass_eV 510998.903           //!< electron mass in eV
+#define Units2Electronmass (AtomicMassUnitsToeV/ElectronMass_eV) //!< atomic mass unit in eV/ electron mass in eV = 1 822.88863
+#define Atomictime2Femtoseconds 0.024188843     //!< Atomictime in fs
+
 #define VERSIONSTRING "v1.0"
 

src/graph.cpp

@@ -7 +7 @@
 using namespace std;
 
+#include "graph.hpp"
 
-#include <iostream>
-#include <list>
-#include <vector>
+/***************************************** Implementations for graph classes ********************************/
 
-/***************************************** Functions for class graph ********************************/
+/** Constructor of class Graph.
+ */
+Graph::Graph()
+{
+};
 
+/** Destructor of class Graph.
+ * Destructor does release memory for nodes and edges contained in its lists as well.
+ */
+Graph::~Graph()
+{
+};
 
+/** Constructor of class SubGraph.
+ */
+SubGraph::SubGraph()
+{
+};
+
+/** Destructor of class SubGraph.
+ * Note that destructor does not deallocate either nodes or edges! (this is done by its subgraph!)
+ */
+SubGraph::~SubGraph()
+{
+};
+
+/** Constructor of class Node.
+ */
+Node::Node()
+{
+};
+
+/** Destructor of class Node.
+ */
+Node::~Node()
+{
+};
+
+/** Constructor of class Edge.
+ */
+Edge::Edge()
+{
+};
+
+/** Destructor of class Edge.
+ */
+Edge::~Edge()
+{
+};
+

src/joiner.cpp

@@ -19 +19 @@
   periodentafel *periode = NULL; // and a period table of all elements
   EnergyMatrix Energy;
+  EnergyMatrix EnergyFragments;
+  
   EnergyMatrix Hcorrection;
+  EnergyMatrix HcorrectionFragments;
+  
   ForceMatrix Force;
-  EnergyMatrix EnergyFragments;
-  EnergyMatrix HcorrectionFragments;
   ForceMatrix ForceFragments;
+
+  HessianMatrix Hessian;
+  HessianMatrix HessianFragments;
+  
   ForceMatrix Shielding;
   ForceMatrix ShieldingPAS;
   ForceMatrix ShieldingFragments;
   ForceMatrix ShieldingPASFragments;
-  KeySetsContainer KeySet;
+  ForceMatrix Chi;
+  ForceMatrix ChiPAS;
+  ForceMatrix ChiFragments;
+  ForceMatrix ChiPASFragments;
+  KeySetsContainer KeySet;  
   stringstream prefix;
   char *dir = NULL;
-  bool Hcorrected = true;
+  bool NoHCorrection = false;
+  bool NoHessian = false;
 
   cout << "Joiner" << endl;
@@ -61 +72 @@
   // ------------- Parse through all Fragment subdirs --------
   if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix, 0,0)) return 1;
-  Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0);
+  if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0)) {
+    NoHCorrection = true;
+    cout << "No HCorrection matrices found, skipping these." << endl;
+  }
   if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1;
+  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix, 0,0)) {
+    NoHessian = true;
+    cout << "No hessian matrices found, skipping these." << endl;
+  }
   if (periode != NULL) { // also look for PAS values
     if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
     if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
+    if (!Chi.ParseFragmentMatrix(argv[1], dir, ChiSuffix, 1, 0)) return 1;
+    if (!ChiPAS.ParseFragmentMatrix(argv[1], dir, ChiPASSuffix, 1, 0)) return 1;
   }
 
   // ---------- Parse the TE Factors into an array -----------------
-  if (!Energy.ParseIndices()) return 1;
-  if (Hcorrected) Hcorrection.ParseIndices();
-
+  if (!Energy.InitialiseIndices()) return 1;
+  if (!NoHCorrection) 
+    Hcorrection.InitialiseIndices();
+  
   // ---------- Parse the Force indices into an array ---------------
   if (!Force.ParseIndices(argv[1])) return 1;
 
+  // ---------- Parse the Hessian (=force) indices into an array ---------------
+  if (!NoHessian)
+    if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
+
   // ---------- Parse the shielding indices into an array ---------------
   if (periode != NULL) { // also look for PAS values
     if(!Shielding.ParseIndices(argv[1])) return 1;
     if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
+    if(!Chi.ParseIndices(argv[1])) return 1;
+    if(!ChiPAS.ParseIndices(argv[1])) return 1;
   }
 
@@ -85 +112 @@
 
   if (!KeySet.ParseManyBodyTerms()) return 1;
+
   if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return 1;
-  if (Hcorrected)  HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter);
+  if (!NoHCorrection)  
+    HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter);
   if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
+  if (!NoHessian)
+    if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1;
   if (periode != NULL) { // also look for PAS values
     if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1;
     if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1;
+    if (!ChiFragments.AllocateMatrix(Chi.Header, Chi.MatrixCounter, Chi.RowCounter, Chi.ColumnCounter)) return 1;
+    if (!ChiPASFragments.AllocateMatrix(ChiPAS.Header, ChiPAS.MatrixCounter, ChiPAS.RowCounter, ChiPAS.ColumnCounter)) return 1;
   }
 
@@ -96 +129 @@
   if(!Energy.SetLastMatrix(0., 0)) return 1;
   if(!Force.SetLastMatrix(0., 2)) return 1;
+  if (!NoHessian)
+    if (!Hessian.SetLastMatrix(0., 0)) return 1;
   if (periode != NULL) { // also look for PAS values
     if(!Shielding.SetLastMatrix(0., 2)) return 1;
     if(!ShieldingPAS.SetLastMatrix(0., 2)) return 1;
+    if(!Chi.SetLastMatrix(0., 2)) return 1;
+    if(!ChiPAS.SetLastMatrix(0., 2)) return 1;
   }
 
@@ -108 +145 @@
     cout << "Summing energy of order " << BondOrder+1 << " ..." << endl;
     if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return 1;
-    if (Hcorrected) {
+    if (!NoHCorrection) { 
       HcorrectionFragments.SumSubManyBodyTerms(Hcorrection, KeySet, BondOrder);
       if (!Energy.SumSubEnergy(EnergyFragments, &HcorrectionFragments, KeySet, BondOrder, 1.)) return 1;
-      if (Hcorrected) Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.);
-    } else
+      Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.);
+    } else 
       if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1;
     // --------- sum up Forces --------------------
@@ -118 +155 @@
     if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return 1;
     if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return 1;
+    // --------- sum up Hessian --------------------
+    if (!NoHessian) {
+      cout << "Summing Hessian of order " << BondOrder+1 << " ..." << endl;
+      if (!HessianFragments.SumSubManyBodyTerms(Hessian, KeySet, BondOrder)) return 1;
+      if (!Hessian.SumSubHessians(HessianFragments, KeySet, BondOrder, 1.)) return 1;
+    }
     if (periode != NULL) { // also look for PAS values
-      cout << "Summing shieldings of order " << BondOrder+1 << " ..." << endl;
+      cout << "Summing shieldings and susceptibilities of order " << BondOrder+1 << " ..." << endl;
       if (!ShieldingFragments.SumSubManyBodyTerms(Shielding, KeySet, BondOrder)) return 1;
       if (!Shielding.SumSubForces(ShieldingFragments, KeySet, BondOrder, 1.)) return 1;
       if (!ShieldingPASFragments.SumSubManyBodyTerms(ShieldingPAS, KeySet, BondOrder)) return 1;
       if (!ShieldingPAS.SumSubForces(ShieldingPASFragments, KeySet, BondOrder, 1.)) return 1;
+      if (!ChiFragments.SumSubManyBodyTerms(Chi, KeySet, BondOrder)) return 1;
+      if (!Chi.SumSubForces(ChiFragments, KeySet, BondOrder, 1.)) return 1;
+      if (!ChiPASFragments.SumSubManyBodyTerms(ChiPAS, KeySet, BondOrder)) return 1;
+      if (!ChiPAS.SumSubForces(ChiPASFragments, KeySet, BondOrder, 1.)) return 1;
     }
 
@@ -134 +181 @@
     // forces
     if (!Force.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ForcesSuffix)) return 1;
+    // hessian
+    if (!NoHessian)
+      if (!Hessian.WriteLastMatrix(argv[1], (prefix.str()).c_str(), HessianSuffix)) return 1;
     // shieldings
     if (periode != NULL) { // also look for PAS values
       if (!Shielding.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ShieldingSuffix)) return 1;
       if (!ShieldingPAS.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ShieldingPASSuffix)) return 1;
+      if (!Chi.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ChiSuffix)) return 1;
+      if (!ChiPAS.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ChiPASSuffix)) return 1;
     }
   }
@@ -144 +196 @@
   prefix << dir << EnergyFragmentSuffix;
   if (!EnergyFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
-  if (Hcorrected) {
+  if (!NoHCorrection) {
     prefix.str(" ");
     prefix << dir << HcorrectionFragmentSuffix;
@@ -153 +205 @@
   if (!ForceFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
   if (!CreateDataFragment(EnergyFragments, KeySet, argv[1], FRAGMENTPREFIX ENERGYPERFRAGMENT, "fragment energy versus the Fragment No", "today", CreateEnergy)) return 1;
+  if (!NoHessian) {
+    prefix.str(" ");
+    prefix << dir << HessianFragmentSuffix;
+    if (!HessianFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+  }
   if (periode != NULL) { // also look for PAS values
     prefix.str(" ");
@@ -160 +217 @@
     prefix << dir << ShieldingPASFragmentSuffix;
     if (!ShieldingPASFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+    prefix.str(" ");
+    prefix << dir << ChiFragmentSuffix;
+    if (!ChiFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+    prefix.str(" ");
+    prefix << dir << ChiPASFragmentSuffix;
+    if (!ChiPASFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
   }
 
   // write last matrices as fragments into central dir (not subdir as above), for analyzer to know index bounds
   if (!Energy.WriteLastMatrix(argv[1], dir, EnergyFragmentSuffix)) return 1;
-  if (Hcorrected) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix);
+  if (!NoHCorrection) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix);
   if (!Force.WriteLastMatrix(argv[1], dir, ForceFragmentSuffix)) return 1;
+  if (!NoHessian)
+    if (!Hessian.WriteLastMatrix(argv[1], dir, HessianFragmentSuffix)) return 1;
   if (periode != NULL) { // also look for PAS values
     if (!Shielding.WriteLastMatrix(argv[1], dir, ShieldingFragmentSuffix)) return 1;
     if (!ShieldingPAS.WriteLastMatrix(argv[1], dir, ShieldingPASFragmentSuffix)) return 1;
+    if (!Chi.WriteLastMatrix(argv[1], dir, ChiFragmentSuffix)) return 1;
+    if (!ChiPAS.WriteLastMatrix(argv[1], dir, ChiPASFragmentSuffix)) return 1;
   }
 

src/moleculelist.cpp

@@ -35 +35 @@
  * \return true - add successful
  */
-bool MoleculeListClass::insert(molecule *mol)
+void MoleculeListClass::insert(molecule *mol)
 {
   mol->IndexNr = MaxIndex++;
   ListOfMolecules.push_back(mol);
-  return true;
 };
 
@@ -128 +127 @@
   atom *Walker = NULL;
   int Counts[MAX_ELEMENTS];
+  double size=0;
+  Vector Origin;
 
   // header
-  *out << "Index\tName\tNo.Atoms\tformula" << endl;
+  *out << "Index\tName\t\tAtoms\tFormula\tCenter\tSize" << endl;
   cout << Verbose(0) << "-----------------------------------------------" << endl;
   if (ListOfMolecules.size() == 0)
     *out << "\tNone" << endl;
   else {
+    Origin.Zero();
     for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
       // reset element counts
       for (int j = 0; j<MAX_ELEMENTS;j++)
         Counts[j] = 0;
-      // count atoms per element
+      // count atoms per element and determine size of bounding sphere
+      size=0.;
       Walker = (*ListRunner)->start;
       while (Walker->next != (*ListRunner)->end) {
         Walker = Walker->next;
         Counts[Walker->type->Z]++;
+        if (Walker->x.DistanceSquared(&Origin) > size)
+          size = Walker->x.DistanceSquared(&Origin);
       }
       // output Index, Name, number of atoms, chemical formula
       *out << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->AtomCount << "\t";
       Elemental = (*ListRunner)->elemente->end;
-      while(Elemental != (*ListRunner)->elemente->start) {
+      while(Elemental->previous != (*ListRunner)->elemente->start) {
         Elemental = Elemental->previous;
         if (Counts[Elemental->Z] != 0)
           *out << Elemental->symbol << Counts[Elemental->Z];
       }
-      *out << endl;
+      // Center and size
+      *out << "\t" << (*ListRunner)->Center << "\t" << sqrt(size) << endl;
     }
   }
@@ -164 +170 @@
 molecule * MoleculeListClass::ReturnIndex(int index)
 {
-  int count = 1;
-  MoleculeList::iterator ListRunner = ListOfMolecules.begin();
-  for(; ((ListRunner != ListOfMolecules.end()) && (count < index)); ListRunner++);
-  if (count == index)
-    return (*ListRunner);
-  else
-    return NULL;
+  for(MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
+    if ((*ListRunner)->IndexNr == index)
+      return (*ListRunner);
+  return NULL;
 };
 
@@ -565 +568 @@
  * \param *configuration standard configuration to attach atoms in fragment molecule to.
  * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
+ * \param DoPeriodic true - call ScanForPeriodicCorrection, false - don't
+ * \param DoCentering true - call molecule::CenterEdge(), false - don't
  * \return true - success (each file was written), false - something went wrong.
  */
-bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out,
-    config *configuration, int *SortIndex)
+bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
 {
   ofstream outputFragment;

src/molecules.cpp

@@ -63 +63 @@
   cell_size[1] = cell_size[3] = cell_size[4]= 0.;
   strcpy(name,"none");
+  IndexNr  = -1;
+  ActiveFlag = false;
 };
 
@@ -602 +604 @@
  * \param *filename filename
  */
-void molecule::SetNameFromFilename(char *filename)
+void molecule::SetNameFromFilename(const char *filename)
 {
   int length = 0;
   char *molname = strrchr(filename, '/')+sizeof(char);  // search for filename without dirs
-  char *endname = strrchr(filename, '.');
+  char *endname = strchr(molname, '.');
   if ((endname == NULL) || (endname < molname))
     length = strlen(molname);
@@ -612 +614 @@
     length = strlen(molname) - strlen(endname);
   strncpy(name, molname, length);
+  name[length]='\0';
 };
 
@@ -669 +672 @@
         }
       }
-      // matrix multiply
-      x.MatrixMultiplication(M);
-      ptr->x.CopyVector(&x);
     }
   }
@@ -710 +710 @@
     max->AddVector(min);
     Translate(min);
+    Center.Zero();
   }
   delete(min);
@@ -719 +720 @@
  * \param *center return vector for translation vector
  */
-void molecule::CenterOrigin(ofstream *out, Vector *center)
+void molecule::CenterOrigin(ofstream *out)
 {
   int Num = 0;
   atom *ptr = start->next;  // start at first in list
 
-  for(int i=NDIM;i--;) // zero center vector
-    center->x[i] = 0.;
+  Center.Zero();
 
   if (ptr != end) {   //list not empty?
@@ -731 +731 @@
       ptr = ptr->next;
       Num++;
-      center->AddVector(&ptr->x);
-    }
-    center->Scale(-1./Num); // divide through total number (and sign for direction)
-    Translate(center);
+      Center.AddVector(&ptr->x);
+    }
+    Center.Scale(-1./Num); // divide through total number (and sign for direction)
+    Translate(&Center);
+    Center.Zero();
   }
 };
@@ -799 +800 @@
  * \param *center return vector for translation vector
  */
-void molecule::CenterGravity(ofstream *out, Vector *center)
-{
-  if (center == NULL) {
-    DetermineCenter(*center);
-    Translate(center);
-    delete(center);
-  } else {
-    Translate(center);
-  }
-};
+void molecule::CenterPeriodic(ofstream *out)
+{
+  DeterminePeriodicCenter(Center);
+};
+
+
+/** Centers the center of gravity of the atoms at (0,0,0).
+ * \param *out output stream for debugging
+ * \param *center return vector for translation vector
+ */
+void molecule::CenterAtVector(ofstream *out, Vector *newcenter)
+{
+  Center.CopyVector(newcenter);
+};
+
 
 /** Scales all atoms by \a *factor.
@@ -911 +917 @@
 
 /** Determines center of molecule (yet not considering atom masses).
- * \param Center reference to return vector
- */
-void molecule::DetermineCenter(Vector &Center)
+ * \param center reference to return vector
+ */
+void molecule::DeterminePeriodicCenter(Vector &center)
 {
   atom *Walker = start;
@@ -987 +993 @@
   Vector *CenterOfGravity = DetermineCenterOfGravity(out);
 
-  CenterGravity(out, CenterOfGravity);
+  CenterPeriodic(out);
 
   // reset inertia tensor
@@ -1083 +1089 @@
 };
 
+/** Evaluates the potential energy used for constrained molecular dynamics.
+ * \f$V_i^{con} = c^{bond} \cdot | r_{P(i)} - R_i | + sum_{i \neq j} C^{min} \cdot \frac{1}{C_{ij}} + C^{inj} \Bigl (1 - \theta \bigl (\prod_{i \neq j} (P(i) - P(j)) \bigr ) \Bigr )\f$
+ *     where the first term points to the target in minimum distance, the second is a penalty for trajectories lying too close to each other (\f$C_{ij}$ is minimum distance between 
+ *     trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
+ * Note that for the second term we have to solve the following linear system:
+ * \f$-c_1 \cdot n_1 + c_2 \cdot n_2 + C \cdot n_3 = - p_2 + p_1\f$, where \f$c_1\f$, \f$c_2\f$ and \f$C\f$ are constants, 
+ * offset vector \f$p_1\f$ in direction \f$n_1\f$, offset vector \f$p_2\f$ in direction \f$n_2\f$,
+ * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines.
+ * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
+ * \param *out output stream for debugging
+ * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
+ * \param startstep start configuration (MDStep in molecule::trajectories)
+ * \param endstep end configuration (MDStep in molecule::trajectories)
+ * \param *constants constant in front of each term
+ * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
+ * \return potential energy
+ * \note This routine is scaling quadratically which is not optimal.
+ * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
+ */
+double molecule::ConstrainedPotential(ofstream *out, atom **PermutationMap, int startstep, int endstep, double *constants, bool IsAngstroem)
+{
+  double result = 0., tmp, Norm1, Norm2;
+  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
+  Vector trajectory1, trajectory2, normal, TestVector;
+  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
+  gsl_vector *x = gsl_vector_alloc(NDIM);
+
+  // go through every atom
+  Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    // first term: distance to target
+    Runner = PermutationMap[Walker->nr];   // find target point
+    tmp = (Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)));
+    tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
+    result += constants[0] * tmp;
+    //*out << Verbose(4) << "Adding " << tmp*constants[0] << "." << endl;
+    
+    // second term: sum of distances to other trajectories
+    Runner = start;
+    while (Runner->next != end) {
+      Runner = Runner->next;
+      if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
+        break;
+      // determine normalized trajectories direction vector (n1, n2)
+      Sprinter = PermutationMap[Walker->nr];   // find first target point
+      trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));
+      trajectory1.SubtractVector(&Trajectories[Walker].R.at(startstep));
+      trajectory1.Normalize();
+      Norm1 = trajectory1.Norm();
+      Sprinter = PermutationMap[Runner->nr];   // find second target point
+      trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));
+      trajectory2.SubtractVector(&Trajectories[Runner].R.at(startstep));
+      trajectory2.Normalize();
+      Norm2 = trajectory1.Norm();
+      // check whether either is zero()
+      if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
+        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
+      } else if (Norm1 < MYEPSILON) {
+        Sprinter = PermutationMap[Walker->nr];   // find first target point
+        trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy first offset
+        trajectory1.SubtractVector(&Trajectories[Runner].R.at(startstep));  // subtract second offset
+        trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything 
+        trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
+        tmp = trajectory1.Norm();  // remaining norm is distance
+      } else if (Norm2 < MYEPSILON) {
+        Sprinter = PermutationMap[Runner->nr];   // find second target point
+        trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy second offset
+        trajectory2.SubtractVector(&Trajectories[Walker].R.at(startstep));  // subtract first offset
+        trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything 
+        trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
+        tmp = trajectory2.Norm();  // remaining norm is distance
+      } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
+//        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
+//        *out << trajectory1;
+//        *out << " and ";
+//        *out << trajectory2;
+        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
+//        *out << " with distance " << tmp << "." << endl;
+      } else { // determine distance by finding minimum distance
+//        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
+//        *out << endl;
+//        *out << "First Trajectory: ";
+//        *out << trajectory1 << endl;
+//        *out << "Second Trajectory: ";
+//        *out << trajectory2 << endl;
+        // determine normal vector for both
+        normal.MakeNormalVector(&trajectory1, &trajectory2);
+        // print all vectors for debugging
+//        *out << "Normal vector in between: ";
+//        *out << normal << endl;
+        // setup matrix
+        for (int i=NDIM;i--;) {
+          gsl_matrix_set(A, 0, i, trajectory1.x[i]);
+          gsl_matrix_set(A, 1, i, trajectory2.x[i]);
+          gsl_matrix_set(A, 2, i, normal.x[i]);
+          gsl_vector_set(x,i, (Trajectories[Walker].R.at(startstep).x[i] - Trajectories[Runner].R.at(startstep).x[i]));
+        }
+        // solve the linear system by Householder transformations
+        gsl_linalg_HH_svx(A, x);
+        // distance from last component
+        tmp = gsl_vector_get(x,2);
+//        *out << " with distance " << tmp << "." << endl;
+        // test whether we really have the intersection (by checking on c_1 and c_2)
+        TestVector.CopyVector(&Trajectories[Runner].R.at(startstep));
+        trajectory2.Scale(gsl_vector_get(x,1));
+        TestVector.AddVector(&trajectory2);
+        normal.Scale(gsl_vector_get(x,2));
+        TestVector.AddVector(&normal);
+        TestVector.SubtractVector(&Trajectories[Walker].R.at(startstep));
+        trajectory1.Scale(gsl_vector_get(x,0));
+        TestVector.SubtractVector(&trajectory1);
+        if (TestVector.Norm() < MYEPSILON) {
+//          *out << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;  
+        } else {
+//          *out << Verbose(2) << "Test: failed.\tIntersection is off by ";
+//          *out << TestVector;
+//          *out << "." << endl;  
+        }
+      }
+      // add up
+      tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
+      if (fabs(tmp) > MYEPSILON) {
+        result += constants[1] * 1./tmp;
+        //*out << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl;
+      }
+    }
+      
+    // third term: penalty for equal targets
+    Runner = start;
+    while (Runner->next != end) {
+      Runner = Runner->next;
+      if ((PermutationMap[Walker->nr] == PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
+        Sprinter = PermutationMap[Walker->nr];
+//        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
+//        *out << Trajectories[Sprinter].R.at(endstep);
+//        *out << ", penalting." << endl;
+        result += constants[2];
+        //*out << Verbose(4) << "Adding " << constants[2] << "." << endl;
+      }
+    }
+  }
+  
+  return result;
+};
+
+void PrintPermutationMap(ofstream *out, atom **PermutationMap, int Nr) 
+{
+  stringstream zeile1, zeile2;
+  int *DoubleList = (int *) Malloc(Nr*sizeof(int), "PrintPermutationMap: *DoubleList");
+  int doubles = 0;
+  for (int i=0;i<Nr;i++)
+    DoubleList[i] = 0;
+  zeile1 << "PermutationMap: ";
+  zeile2 << "                ";
+  for (int i=0;i<Nr;i++) {
+    DoubleList[PermutationMap[i]->nr]++;
+    zeile1 << i << " ";
+    zeile2 << PermutationMap[i]->nr << " ";
+  }
+  for (int i=0;i<Nr;i++)
+    if (DoubleList[i] > 1)
+    doubles++;
+ // *out << "Found " << doubles << " Doubles." << endl;
+  Free((void **)&DoubleList, "PrintPermutationMap: *DoubleList");
+//  *out << zeile1.str() << endl << zeile2.str() << endl;
+};
+
+/** Minimises the extra potential for constrained molecular dynamics and gives forces and the constrained potential energy.
+ * We do the following:
+ *  -# Generate a distance list from all source to all target points
+ *  -# Sort this per source point
+ *  -# Take for each source point the target point with minimum distance, use this as initial permutation
+ *  -# check whether molecule::ConstrainedPotential() is greater than injective penalty
+ *     -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential.
+ *  -# Next, we only apply transformations that keep the injectivity of the permutations list.
+ *  -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target 
+ *     point and try to change it for one with lesser distance, or for the next one with greater distance, but only
+ *     if this decreases the conditional potential.
+ *  -# finished. 
+ *  -# Then, we calculate the forces by taking the spatial derivative, where we scale the potential to such a degree,
+ *     that the total force is always pointing in direction of the constraint force (ensuring that we move in the
+ *     right direction).
+ *  -# Finally, we calculate the potential energy and return.
+ * \param *out output stream for debugging
+ * \param **PermutationMap on return: mapping between the atom label of the initial and the final configuration
+ * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
+ * \param endstep step giving final position in constrained MD 
+ * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
+ * \sa molecule::VerletForceIntegration()
+ * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2)
+ * \todo The constrained potential's constants are set to fixed values right now, but they should scale based on checks of the system in order
+ *       to ensure they're properties (e.g. constants[2] always greater than the energy of the system).
+ * \bug this all is not O(N log N) but O(N^2)
+ */
+double molecule::MinimiseConstrainedPotential(ofstream *out, atom **&PermutationMap, int startstep, int endstep, bool IsAngstroem)
+{
+  double Potential, OldPotential, OlderPotential;
+  PermutationMap = (atom **) Malloc(AtomCount*sizeof(atom *), "molecule::MinimiseConstrainedPotential: **PermutationMap");
+  DistanceMap **DistanceList = (DistanceMap **) Malloc(AtomCount*sizeof(DistanceMap *), "molecule::MinimiseConstrainedPotential: **DistanceList");
+  DistanceMap::iterator *DistanceIterators = (DistanceMap::iterator *) Malloc(AtomCount*sizeof(DistanceMap::iterator), "molecule::MinimiseConstrainedPotential: *DistanceIterators");
+  int *DoubleList = (int *) Malloc(AtomCount*sizeof(int), "molecule::MinimiseConstrainedPotential: *DoubleList");
+  DistanceMap::iterator *StepList = (DistanceMap::iterator *) Malloc(AtomCount*sizeof(DistanceMap::iterator), "molecule::MinimiseConstrainedPotential: *StepList");
+  double constants[3];
+  int round;
+  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
+  DistanceMap::iterator Rider, Strider;
+  
+  /// Minimise the potential
+  // set Lagrange multiplier constants
+  constants[0] = 10.;
+  constants[1] = 1.;
+  constants[2] = 1e+7;    // just a huge penalty
+  // generate the distance list
+  *out << Verbose(1) << "Creating the distance list ... " << endl;
+  for (int i=AtomCount; i--;) {
+    DoubleList[i] = 0;  // stores for how many atoms in startstep this atom is a target in endstep
+    DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
+    DistanceList[i]->clear();
+  }
+  *out << Verbose(1) << "Filling the distance list ... " << endl;
+  Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    Runner = start;
+    while(Runner->next != end) {
+      Runner = Runner->next;
+      DistanceList[Walker->nr]->insert( DistancePair(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)), Runner) );
+    }
+  }
+  // create the initial PermutationMap (source -> target)
+  Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    StepList[Walker->nr] = DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
+    PermutationMap[Walker->nr] = DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
+    DoubleList[DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
+    DistanceIterators[Walker->nr] = DistanceList[Walker->nr]->begin();    // and remember which one we picked
+    *out << *Walker << " starts with distance " << DistanceList[Walker->nr]->begin()->first << "." << endl;
+  }
+  *out << Verbose(1) << "done." << endl;
+  // make the PermutationMap injective by checking whether we have a non-zero constants[2] term in it
+  *out << Verbose(1) << "Making the PermutationMap injective ... " << endl;
+  Walker = start;
+  DistanceMap::iterator NewBase;
+  OldPotential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
+  while ((OldPotential) > constants[2]) {
+    PrintPermutationMap(out, PermutationMap, AtomCount);
+    Walker = Walker->next;
+    if (Walker == end) // round-robin at the end
+      Walker = start->next;
+    if (DoubleList[DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
+      continue;
+    // now, try finding a new one
+    NewBase = DistanceIterators[Walker->nr];  // store old base
+    do {
+      NewBase++;  // take next further distance in distance to targets list that's a target of no one
+    } while ((DoubleList[NewBase->second->nr] != 0) && (NewBase != DistanceList[Walker->nr]->end()));
+    if (NewBase != DistanceList[Walker->nr]->end()) {
+      PermutationMap[Walker->nr] = NewBase->second;
+      Potential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
+      if (Potential > OldPotential) { // undo
+        PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
+      } else {  // do
+        DoubleList[DistanceIterators[Walker->nr]->second->nr]--;  // decrease the old entry in the doubles list
+        DoubleList[NewBase->second->nr]++;    // increase the old entry in the doubles list
+        DistanceIterators[Walker->nr] = NewBase;
+        OldPotential = Potential;
+        *out << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl;
+      }
+    }
+  } 
+  for (int i=AtomCount; i--;) // now each single entry in the DoubleList should be <=1
+    if (DoubleList[i] > 1) { 
+      cerr << "Failed to create an injective PermutationMap!" << endl;
+      exit(1);
+    }
+  *out << Verbose(1) << "done." << endl;
+  Free((void **)&DoubleList, "molecule::MinimiseConstrainedPotential: *DoubleList");
+  // argument minimise the constrained potential in this injective PermutationMap
+  *out << Verbose(1) << "Argument minimising the PermutationMap, at current potential " << OldPotential << " ... " << endl;
+  OldPotential = 1e+10;
+  round = 0;
+  do {
+    *out << "Starting round " << ++round << " ... " << endl;
+    OlderPotential = OldPotential;
+    do {
+      Walker = start;
+      while (Walker->next != end) { // pick one
+        Walker = Walker->next;
+        PrintPermutationMap(out, PermutationMap, AtomCount);
+        Sprinter = DistanceIterators[Walker->nr]->second;   // store initial partner
+        Strider = DistanceIterators[Walker->nr];  //remember old iterator
+        DistanceIterators[Walker->nr] = StepList[Walker->nr];  
+        if (DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
+          DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->begin();
+          break;
+        }
+        //*out << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl;
+        // find source of the new target
+        Runner = start->next;
+        while(Runner != end) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
+          if (PermutationMap[Runner->nr] == DistanceIterators[Walker->nr]->second) {
+            //*out << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *PermutationMap[Runner->nr] << "." << endl;
+            break;
+          }
+          Runner = Runner->next;
+        }
+        if (Runner != end) { // we found the other source
+          // then look in its distance list for Sprinter
+          Rider = DistanceList[Runner->nr]->begin();
+          for (; Rider != DistanceList[Runner->nr]->end(); Rider++)
+            if (Rider->second == Sprinter)
+              break;
+          if (Rider != DistanceList[Runner->nr]->end()) { // if we have found one
+            //*out << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl; 
+            // exchange both 
+            PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap 
+            PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
+            PrintPermutationMap(out, PermutationMap, AtomCount);
+            // calculate the new potential
+            //*out << Verbose(2) << "Checking new potential ..." << endl;
+            Potential = ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
+            if (Potential > OldPotential) { // we made everything worse! Undo ...
+              //*out << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
+              //*out << Verbose(3) << "Setting " << *Runner << "'s source to " << *DistanceIterators[Runner->nr]->second << "." << endl; 
+              // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
+              PermutationMap[Runner->nr] = DistanceIterators[Runner->nr]->second;
+              // Undo for Walker
+              DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
+              //*out << Verbose(3) << "Setting " << *Walker << "'s source to " << *DistanceIterators[Walker->nr]->second << "." << endl; 
+              PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
+            } else {
+              DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
+              *out << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl;
+              OldPotential = Potential;
+            }
+            if (Potential > constants[2]) {
+              cerr << "ERROR: The two-step permutation procedure did not maintain injectivity!" << endl;
+              exit(255);
+            }
+            //*out << endl;
+          } else {
+            cerr << "ERROR: " << *Runner << " was not the owner of " << *Sprinter << "!" << endl;
+            exit(255);
+          }
+        } else {
+          PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // new target has no source!
+        }
+        StepList[Walker->nr]++; // take next farther distance target
+      }
+    } while (Walker->next != end);
+  } while ((OlderPotential - OldPotential) > 1e-3);
+  *out << Verbose(1) << "done." << endl;
+
+  
+  /// free memory and return with evaluated potential 
+  for (int i=AtomCount; i--;)
+    DistanceList[i]->clear();
+  Free((void **)&DistanceList, "molecule::MinimiseConstrainedPotential: **DistanceList");
+  Free((void **)&DistanceIterators, "molecule::MinimiseConstrainedPotential: *DistanceIterators");
+  return ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
+};
+
+/** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
+ * \param *out output stream for debugging
+ * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
+ * \param endstep step giving final position in constrained MD 
+ * \param **PermutationMap mapping between the atom label of the initial and the final configuration
+ * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
+ * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
+ */
+void molecule::EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
+{
+  double constant = 10.;
+  atom *Walker = NULL, *Sprinter = NULL; 
+  
+  /// evaluate forces (only the distance to target dependent part) with the final PermutationMap
+  *out << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl;
+  Walker = start;
+  while (Walker->next != NULL) {
+    Walker = Walker->next;
+    Sprinter = PermutationMap[Walker->nr];
+    // set forces
+    for (int i=NDIM;i++;)
+      Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Sprinter].R.at(endstep)));
+  }  
+  *out << Verbose(1) << "done." << endl;
+};
+
+/** Performs a linear interpolation between two desired atomic configurations with a given number of steps.
+ * Note, step number is config::MaxOuterStep
+ * \param *out output stream for debugging
+ * \param startstep stating initial configuration in molecule::Trajectories
+ * \param endstep stating final configuration in molecule::Trajectories
+ * \param &config configuration structure
+ * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
+ */
+bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration) 
+{
+  molecule *mol = NULL;
+  bool status = true;
+  int MaxSteps = configuration.MaxOuterStep; 
+  MoleculeListClass *MoleculePerStep = new MoleculeListClass();
+  // Get the Permutation Map by MinimiseConstrainedPotential
+  atom **PermutationMap = NULL;
+  atom *Walker = NULL, *Sprinter = NULL;
+  MinimiseConstrainedPotential(out, PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
+
+  // check whether we have sufficient space in Trajectories for each atom
+  Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if (Trajectories[Walker].R.size() <= (unsigned int)(MaxSteps)) {
+      //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
+      Trajectories[Walker].R.resize(MaxSteps+1);
+      Trajectories[Walker].U.resize(MaxSteps+1);
+      Trajectories[Walker].F.resize(MaxSteps+1);
+    }
+  }
+  // push endstep to last one
+  Walker = start;
+  while (Walker->next != end) { // remove the endstep (is now the very last one)
+    Walker = Walker->next;
+    for (int n=NDIM;n--;) {
+      Trajectories[Walker].R.at(MaxSteps).x[n] = Trajectories[Walker].R.at(endstep).x[n];
+      Trajectories[Walker].U.at(MaxSteps).x[n] = Trajectories[Walker].U.at(endstep).x[n];
+      Trajectories[Walker].F.at(MaxSteps).x[n] = Trajectories[Walker].F.at(endstep).x[n];
+    }
+  }
+  endstep = MaxSteps;
+  
+  // go through all steps and add the molecular configuration to the list and to the Trajectories of \a this molecule
+  *out << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl;
+  for (int step = 0; step <= MaxSteps; step++) {
+    mol = new molecule(elemente);
+    MoleculePerStep->insert(mol);
+    Walker = start;
+    while (Walker->next != end) {
+      Walker = Walker->next;
+      // add to molecule list
+      Sprinter = mol->AddCopyAtom(Walker);
+      for (int n=NDIM;n--;) {
+        Sprinter->x.x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
+        // add to Trajectories
+        //*out << Verbose(3) << step << ">=" << MDSteps-1 << endl;
+        if (step < MaxSteps) {
+          Trajectories[Walker].R.at(step).x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
+          Trajectories[Walker].U.at(step).x[n] = 0.;
+          Trajectories[Walker].F.at(step).x[n] = 0.;
+        }
+      }
+    }
+  }
+  MDSteps = MaxSteps+1;   // otherwise new Trajectories' points aren't stored on save&exit
+  
+  // store the list to single step files
+  int *SortIndex = (int *) Malloc(AtomCount*sizeof(int), "molecule::LinearInterpolationBetweenConfiguration: *SortIndex");
+  for (int i=AtomCount; i--; )
+    SortIndex[i] = i;
+  status = MoleculePerStep->OutputConfigForListOfFragments(out, &configuration, SortIndex);
+  
+  // free and return
+  Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
+  delete(MoleculePerStep);
+  return status;
+};
+
 /** Parses nuclear forces from file and performs Verlet integration.
  * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
  * have to transform them).
  * This adds a new MD step to the config file.
+ * \param *out output stream for debugging
  * \param *file filename
+ * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
  * \param delta_t time step width in atomic units
  * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
+ * \param DoConstrained whether we perform a constrained (>0, target step in molecule::trajectories) or unconstrained (0) molecular dynamics, \sa molecule::MinimiseConstrainedPotential() 
  * \return true - file found and parsed, false - file not found or imparsable
- */
-bool molecule::VerletForceIntegration(char *file, double delta_t, bool IsAngstroem)
-{
-  element *runner = elemente->start;
+ * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
+ */
+bool molecule::VerletForceIntegration(ofstream *out, char *file, config &configuration)
+{
   atom *walker = NULL;
-  int AtomNo;
   ifstream input(file);
   string token;
   stringstream item;
-  double a, IonMass;
+  double IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp;
   ForceMatrix Force;
-  Vector tmpvector;
 
   CountElements();  // make sure ElementsInMolecule is up to date
@@ -1120 +1595 @@
     }
     // correct Forces
-//    for(int d=0;d<NDIM;d++)
-//      tmpvector.x[d] = 0.;
-//    for(int i=0;i<AtomCount;i++)
-//      for(int d=0;d<NDIM;d++) {
-//        tmpvector.x[d] += Force.Matrix[0][i][d+5];
-//      }
-//    for(int i=0;i<AtomCount;i++)
-//      for(int d=0;d<NDIM;d++) {
-//        Force.Matrix[0][i][d+5] -= tmpvector.x[d]/(double)AtomCount;
-//      }
+    for(int d=0;d<NDIM;d++)
+      Vector[d] = 0.;
+    for(int i=0;i<AtomCount;i++)
+      for(int d=0;d<NDIM;d++) {
+        Vector[d] += Force.Matrix[0][i][d+5];
+      }
+    for(int i=0;i<AtomCount;i++)
+      for(int d=0;d<NDIM;d++) {
+        Force.Matrix[0][i][d+5] -= Vector[d]/(double)AtomCount;
+      }
+    // solve a constrained potential if we are meant to
+    if (configuration.DoConstrainedMD) {
+      // calculate forces and potential
+      atom **PermutationMap = NULL;
+      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
+      EvaluateConstrainedForces(out, configuration.DoConstrainedMD, 0, PermutationMap, &Force);
+      Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
+    }
+    
     // and perform Verlet integration for each atom with position, velocity and force vector
-    runner = elemente->start;
-    while (runner->next != elemente->end) { // go through every element
-      runner = runner->next;
-      IonMass = runner->mass;
-      a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
-      if (ElementsInMolecule[runner->Z]) { // if this element got atoms
-        AtomNo = 0;
+    walker = start;
+    while (walker->next != end) { // go through every atom of this element
+      walker = walker->next;
+      //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
+      // check size of vectors
+      if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
+        //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
+        Trajectories[walker].R.resize(MDSteps+10);
+        Trajectories[walker].U.resize(MDSteps+10);
+        Trajectories[walker].F.resize(MDSteps+10);
+      }
+      
+      // Update R (and F)
+      for (int d=0; d<NDIM; d++) {
+        Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+        Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
+        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2 
+        Trajectories[walker].R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(Trajectories[walker].F.at(MDSteps).x[d]/walker->type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
+      }
+      // Update U
+      for (int d=0; d<NDIM; d++) {
+        Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d]; 
+        Trajectories[walker].U.at(MDSteps).x[d] += configuration.Deltat * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
+      }
+//      out << "Integrated position&velocity of step " << (MDSteps) << ": ("; 
+//      for (int d=0;d<NDIM;d++)
+//        out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
+//      out << ")\t(";
+//      for (int d=0;d<NDIM;d++)
+//        cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
+//      out << ")" << endl;
+            // next atom
+    }
+  }
+  // correct velocities (rather momenta) so that center of mass remains motionless
+  for(int d=0;d<NDIM;d++)
+    Vector[d] = 0.;
+  IonMass = 0.;
+  walker = start;
+  while (walker->next != end) { // go through every atom
+    walker = walker->next;
+    IonMass += walker->type->mass;  // sum up total mass
+    for(int d=0;d<NDIM;d++) {
+      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
+    }
+  }
+  // correct velocities (rather momenta) so that center of mass remains motionless
+  for(int d=0;d<NDIM;d++)
+    Vector[d] /= IonMass;
+  ActualTemp = 0.;
+  walker = start;
+  while (walker->next != end) { // go through every atom of this element
+    walker = walker->next;
+    for(int d=0;d<NDIM;d++) {
+      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d];
+      ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
+    }
+  }
+  Thermostats(configuration, ActualTemp, Berendsen);
+  MDSteps++;
+
+
+  // exit
+  return true;
+};
+
+/** Implementation of various thermostats.
+ * All these thermostats apply an additional force which has the following forms:
+ * -# Woodcock
+ *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
+ * -# Gaussian
+ *  \f$ \frac{ \sum_i \frac{p_i}{m_i} \frac{\partial V}{\partial q_i}} {\sum_i \frac{p^2_i}{m_i}} \cdot p_i\f$
+ * -# Langevin
+ *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$  
+ * -# Berendsen
+ *  \f$p_i \rightarrow \left [ 1+ \frac{\delta t}{\tau_T} \left ( \frac{T_0}{T} \right ) \right ]^{\frac{1}{2}} \cdot p_i\f$
+ * -# Nose-Hoover
+ *  \f$\zeta p_i \f$ with \f$\frac{\partial \zeta}{\partial t} = \frac{1}{M_s} \left ( \sum^N_{i=1} \frac{p_i^2}{m_i} - g k_B T \right )\f$
+ * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
+ * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
+ * belatedly and the constraint force set.
+ * \param *P Problem at hand
+ * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
+ * \sa InitThermostat()
+ */
+void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
+{
+  double ekin = 0.;
+  double E = 0., G = 0.;
+  double delta_alpha = 0.;
+  double ScaleTempFactor;
+  double sigma;
+  double IonMass;
+  int d;
+  gsl_rng * r;
+  const gsl_rng_type * T;
+  double *U = NULL, *F = NULL, FConstraint[NDIM];
+  atom *walker = NULL;
+  
+  // calculate scale configuration
+  ScaleTempFactor = configuration.TargetTemp/ActualTemp;
+    
+  // differentating between the various thermostats
+  switch(Thermostat) {
+     case None:
+      cout << Verbose(2) <<  "Applying no thermostat..." << endl;
+      break;
+     case Woodcock:
+      if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
+        cout << Verbose(2) <<  "Applying Woodcock thermostat..." << endl;
         walker = start;
         while (walker->next != end) { // go through every atom of this element
           walker = walker->next;
-          if (walker->type == runner) { // if this atom fits to element
-            // check size of vectors
-            if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
-              //cout << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
-              Trajectories[walker].R.resize(MDSteps+10);
-              Trajectories[walker].U.resize(MDSteps+10);
-              Trajectories[walker].F.resize(MDSteps+10);
+          IonMass = walker->type->mass;
+          U = Trajectories[walker].U.at(MDSteps).x;
+          if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
+            for (d=0; d<NDIM; d++) {
+              U[d] *= sqrt(ScaleTempFactor);
+              ekin += 0.5*IonMass * U[d]*U[d];
             }
-            // 1. calculate x(t+\delta t)
-            for (int d=0; d<NDIM; d++) {
-              Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5];
-              Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
-              Trajectories[walker].R.at(MDSteps).x[d] += delta_t*(Trajectories[walker].U.at(MDSteps-1).x[d]);
-              Trajectories[walker].R.at(MDSteps).x[d] += 0.5*delta_t*delta_t*(Trajectories[walker].F.at(MDSteps-1).x[d])/IonMass;    // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
+        }
+      }
+      break;
+     case Gaussian:
+      cout << Verbose(2) <<  "Applying Gaussian thermostat..." << endl;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            G += U[d] * F[d]; 
+            E += U[d]*U[d]*IonMass;
+          }
+      }
+      cout << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            FConstraint[d] = (G/E) * (U[d]*IonMass);
+            U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
+            ekin += IonMass * U[d]*U[d];
+          }
+      }
+      break;
+     case Langevin:
+      cout << Verbose(2) <<  "Applying Langevin thermostat..." << endl;
+      // init random number generator
+      gsl_rng_env_setup();
+      T = gsl_rng_default;
+      r = gsl_rng_alloc (T);
+      // Go through each ion
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          // throw a dice to determine whether it gets hit by a heat bath particle
+          if (((((rand()/(double)RAND_MAX))*configuration.TempFrequency) < 1.)) {  
+            cout << Verbose(3) << "Particle " << *walker << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
+            // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
+            for (d=0; d<NDIM; d++) {
+              U[d] = gsl_ran_gaussian (r, sigma);
             }
-            // 2. Calculate v(t+\delta t)
-            for (int d=0; d<NDIM; d++) {
-              Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
-              Trajectories[walker].U.at(MDSteps).x[d] += 0.5*delta_t*(Trajectories[walker].F.at(MDSteps-1).x[d]+Trajectories[walker].F.at(MDSteps).x[d])/IonMass;
-            }
-//            cout << "Integrated position&velocity of step " << (MDSteps) << ": (";
-//            for (int d=0;d<NDIM;d++)
-//              cout << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
-//            cout << ")\t(";
-//            for (int d=0;d<NDIM;d++)
-//              cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
-//            cout << ")" << endl;
-            // next atom
-            AtomNo++;
+            cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
+          }
+          for (d=0; d<NDIM; d++)
+            ekin += 0.5*IonMass * U[d]*U[d];
+        }
+      }
+      break;
+     case Berendsen:
+      cout << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        F = Trajectories[walker].F.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            U[d] *= sqrt(1+(configuration.Deltat/configuration.TempFrequency)*(ScaleTempFactor-1));
+            ekin += 0.5*IonMass * U[d]*U[d];
           }
         }
       }
-    }
-  }
-//  // correct velocities (rather momenta) so that center of mass remains motionless
-//  tmpvector.zero()
-//  IonMass = 0.;
-//  walker = start;
-//  while (walker->next != end) { // go through every atom
-//    walker = walker->next;
-//    IonMass += walker->type->mass;  // sum up total mass
-//    for(int d=0;d<NDIM;d++) {
-//      tmpvector.x[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
-//    }
-//  }
-//  walker = start;
-//  while (walker->next != end) { // go through every atom of this element
-//    walker = walker->next;
-//    for(int d=0;d<NDIM;d++) {
-//      Trajectories[walker].U.at(MDSteps).x[d] -= tmpvector.x[d]*walker->type->mass/IonMass;
-//    }
-//  }
-  MDSteps++;
-
-
-  // exit
-  return true;
-};
-
-/** Align all atoms in such a manner that given vector \a *n is along z axis.
+      break;
+     case NoseHoover:
+      cout << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl;
+      // dynamically evolve alpha (the additional degree of freedom)
+      delta_alpha = 0.;
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+            delta_alpha += U[d]*U[d]*IonMass;
+          }
+        }
+      }
+      delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
+      configuration.alpha += delta_alpha*configuration.Deltat;
+      cout << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl;
+      // apply updated alpha as additional force
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        IonMass = walker->type->mass;
+        U = Trajectories[walker].U.at(MDSteps).x;
+        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
+          for (d=0; d<NDIM; d++) {
+              FConstraint[d] = - configuration.alpha * (U[d] * IonMass);
+              U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
+              ekin += (0.5*IonMass) * U[d]*U[d];
+            }
+        }
+      }
+      break;
+  }    
+  cout << Verbose(1) << "Kinetic energy is " << ekin << "." << endl;
+};
+
+/** Align all atoms in such a manner that given vector \a *n is along z axis. 
  * \param n[] alignment vector.
  */
@@ -1267 +1916 @@
 bool molecule::RemoveAtom(atom *pointer)
 {
-  if (ElementsInMolecule[pointer->type->Z] != 0)  // this would indicate an error
+  if (ElementsInMolecule[pointer->type->Z] != 0)  { // this would indicate an error
     ElementsInMolecule[pointer->type->Z]--;  // decrease number of atom of this element
-  else
+    AtomCount--;
+  } else
     cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
@@ -3093 +3743 @@
   while (MolecularWalker->next != NULL) {
     MolecularWalker = MolecularWalker->next;
+    *out << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
     LocalBackEdgeStack = new StackClass<bond *> (MolecularWalker->Leaf->BondCount);
 //    // check the list of local atoms for debugging
@@ -3108 +3759 @@
     delete(LocalBackEdgeStack);
   }
-
+  
   // ===== 3. if structure still valid, parse key set file and others =====
   FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
@@ -3114 +3765 @@
   // ===== 4. check globally whether there's something to do actually (first adaptivity check)
   FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);
-
-  // =================================== Begin of FRAGMENTATION ===============================
-  // ===== 6a. assign each keyset to its respective subgraph =====
+  
+  // =================================== Begin of FRAGMENTATION =============================== 
+  // ===== 6a. assign each keyset to its respective subgraph ===== 
   Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), true);
 
@@ -3151 +3802 @@
   delete(ParsedFragmentList);
   delete[](MinimumRingSize);
-
+  
 
   // ==================================== End of FRAGMENTATION ============================================
@@ -3244 +3895 @@
   atom *Walker = NULL, *OtherAtom = NULL;
   ReferenceStack->Push(Binder);
-
+  
   do {  // go through all bonds and push local ones
     Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
@@ -3714 +4365 @@
 };
 
+/** Creates \a MoleculeListClass of all unique fragments of the \a molecule containing \a Order atoms or vertices.
+ * The picture to have in mind is that of a DFS "snake" of a certain length \a Order, i.e. as in the infamous
+ * computer game, that winds through the connected graph representing the molecule. Color (white,
+ * lightgray, darkgray, black) indicates whether a vertex has been discovered so far or not. Labels will help in
+ * creating only unique fragments and not additional ones with vertices simply in different sequence.
+ * The Predecessor is always the one that came before in discovering, needed on backstepping. And
+ * finally, the ShortestPath is needed for removing vertices from the snake stack during the back-
+ * stepping.
+ * \param *out output stream for debugging
+ * \param Order number of atoms in each fragment
+ * \param *configuration configuration for writing config files for each fragment
+ * \return List of all unique fragments with \a Order atoms
+ */
+/*
+MoleculeListClass * molecule::CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, config *configuration)
+{
+  atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
+  int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
+  int *Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
+  enum Shading *ColorVertexList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
+  enum Shading *ColorEdgeList = (enum Shading *) Malloc(sizeof(enum Shading)*BondCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorBondList");
+  StackClass<atom *> *RootStack = new StackClass<atom *>(AtomCount);
+  StackClass<atom *> *TouchedStack = new StackClass<atom *>((int)pow(4,Order)+2); // number of atoms reached from one with maximal 4 bonds plus Root itself
+  StackClass<atom *> *SnakeStack = new StackClass<atom *>(Order+1); // equal to Order is not possible, as then the StackClass<atom *> cannot discern between full and empty stack!
+  MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL;
+  MoleculeListClass *FragmentList = NULL;
+  atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
+  bond *Binder = NULL;
+  int RunningIndex = 0, FragmentCounter = 0;
+
+  *out << Verbose(1) << "Begin of CreateListOfUniqueFragmentsOfOrder." << endl;
+
+  // reset parent list
+  *out << Verbose(3) << "Resetting labels, parent, predecessor, color and shortest path lists." << endl;
+  for (int i=0;i<AtomCount;i++) { // reset all atom labels
+    // initialise each vertex as white with no predecessor, empty queue, color lightgray, not labelled, no sons
+    Labels[i] = -1;
+    SonList[i] = NULL;
+    PredecessorList[i] = NULL;
+    ColorVertexList[i] = white;
+    ShortestPathList[i] = -1;
+  }
+  for (int i=0;i<BondCount;i++)
+    ColorEdgeList[i] = white;
+  RootStack->ClearStack();  // clearstack and push first atom if exists
+  TouchedStack->ClearStack();
+  Walker = start->next;
+  while ((Walker != end)
+#ifdef ADDHYDROGEN
+   && (Walker->type->Z == 1)
+        }
+      }
+      *out << ", SP of " << ShortestPathList[Walker->nr]  << " and its color is " << GetColor(ColorVertexList[Walker->nr]) << "." << endl;
+
+      // then check the stack for a newly stumbled upon fragment
+      if (SnakeStack->ItemCount() == Order) { // is stack full?
+        // store the fragment if it is one and get a removal candidate
+        Removal = StoreFragmentFromStack(out, Root, Walker, Leaflet, SnakeStack, ShortestPathList, SonList, Labels, &FragmentCounter, configuration);
+        // remove the candidate if one was found
+        if (Removal != NULL) {
+          *out << Verbose(2) << "Removing item " << Removal->Name << " with SP of " << ShortestPathList[Removal->nr] << " from snake stack." << endl;
+          SnakeStack->RemoveItem(Removal);
+          ColorVertexList[Removal->nr] = lightgray; // return back to not on snake stack but explored marking
+          if (Walker == Removal) { // if the current atom is to be removed, we also have to take a step back
+            Walker = PredecessorList[Removal->nr];
+            *out << Verbose(2) << "Stepping back to " << Walker->Name << "." << endl;
+          }
+        }
+      } else
+        Removal = NULL;
+
+      // finally, look for a white neighbour as the next Walker
+      Binder = NULL;
+      if ((Removal == NULL) || (Walker != PredecessorList[Removal->nr])) {  // don't look, if a new walker has been set above
+        *out << Verbose(2) << "Snake has currently " << SnakeStack->ItemCount() << " item(s)." << endl;
+        OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
+        if (ShortestPathList[Walker->nr] < Order) {
+          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
+            Binder = ListOfBondsPerAtom[Walker->nr][i];
+            *out << Verbose(2) << "Current bond is " << *Binder << ": ";
+            OtherAtom = Binder->GetOtherAtom(Walker);
+            if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
+              *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
+              //ColorVertexList[OtherAtom->nr] = lightgray;    // mark as explored
+            } else { // otherwise check its colour and element
+              if (
+              (OtherAtom->type->Z != 1) &&
+#endif
+                    (ColorEdgeList[Binder->nr] == white)) {  // skip hydrogen, look for unexplored vertices
+                *out << "Moving along " << GetColor(ColorEdgeList[Binder->nr]) << " bond " << Binder << " to " << ((ColorVertexList[OtherAtom->nr] == white) ? "unexplored" : "explored") << " item: " << OtherAtom->Name << "." << endl;
+                // i find it currently rather sensible to always set the predecessor in order to find one's way back
+                //if (PredecessorList[OtherAtom->nr] == NULL) {
+                PredecessorList[OtherAtom->nr] = Walker;
+                *out << Verbose(3) << "Setting Predecessor of " << OtherAtom->Name << " to " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
+                //} else {
+                //  *out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
+                //}
+                Walker = OtherAtom;
+                break;
+              } else {
+                if (OtherAtom->type->Z == 1)
+                  *out << "Links to a hydrogen atom." << endl;
+                else
+                  *out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
+              }
+            }
+          }
+        } else {  // means we have stepped beyond the horizon: Return!
+          Walker = PredecessorList[Walker->nr];
+          OtherAtom = Walker;
+          *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
+        }
+        if (Walker != OtherAtom) {  // if no white neighbours anymore, color it black
+          *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
+          ColorVertexList[Walker->nr] = black;
+          Walker = PredecessorList[Walker->nr];
+        }
+      }
+    } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
+    *out << Verbose(2) << "Inner Looping is finished." << endl;
+
+    // if we reset all AtomCount atoms, we have again technically O(N^2) ...
+    *out << Verbose(2) << "Resetting lists." << endl;
+    Walker = NULL;
+    Binder = NULL;
+    while (!TouchedStack->IsEmpty()) {
+      Walker = TouchedStack->PopLast();
+      *out << Verbose(3) << "Re-initialising entries of " << *Walker << "." << endl;
+      for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
+        ColorEdgeList[ListOfBondsPerAtom[Walker->nr][i]->nr] = white;
+      PredecessorList[Walker->nr] = NULL;
+      ColorVertexList[Walker->nr] = white;
+      ShortestPathList[Walker->nr] = -1;
+    }
+  }
+  *out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
+
+  // copy together
+  *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
+  FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
+  RunningIndex = 0;
+  while ((Leaflet != NULL) && (RunningIndex < FragmentCounter))  {
+    FragmentList->ListOfMolecules[RunningIndex++] = Leaflet->Leaf;
+    Leaflet->Leaf = NULL; // prevent molecule from being removed
+    TempLeaf = Leaflet;
+    Leaflet = Leaflet->previous;
+    delete(TempLeaf);
+  };
+
+  // free memory and exit
+  Free((void **)&PredecessorList, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
+  Free((void **)&ShortestPathList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
+  Free((void **)&Labels, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
+  Free((void **)&ColorVertexList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
+  delete(RootStack);
+  delete(TouchedStack);
+  delete(SnakeStack);
+
+  *out << Verbose(1) << "End of CreateListOfUniqueFragmentsOfOrder." << endl;
+  return FragmentList;
+};
+*/
+
 /** Structure containing all values in power set combination generation.
  */
@@ -4545 +5359 @@
   if (result) {
     *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
-    DetermineCenter(CenterOfGravity);
-    OtherMolecule->DetermineCenter(OtherCenterOfGravity);
+    DeterminePeriodicCenter(CenterOfGravity);
+    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
     *out << Verbose(5) << "Center of Gravity: ";
     CenterOfGravity.Output(out);

src/molecules.hpp

@@ -10 +10 @@
 
 // GSL headers
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_heapsort.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
 #include <gsl/gsl_multimin.h>
 #include <gsl/gsl_vector.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_eigen.h>
-#include <gsl/gsl_heapsort.h>
+#include <gsl/gsl_randist.h>
 
 // STL headers
@@ -156 +158 @@
 };
 
-ostream & operator << (ostream &ost, atom &a);
+ostream & operator << (ostream &ost, const atom &a);
 
 /** Bonds between atoms.
@@ -195 +197 @@
 };
 
-ostream & operator << (ostream &ost, bond &b);
+
+ostream & operator << (ostream &ost, const bond &b);
 
 class MoleculeLeafClass;
+
+
+#define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented 
+enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output
+
 
 /** The complete molecule.
@@ -254 +262 @@
   bool CenterInBox(ofstream *out);
   void CenterEdge(ofstream *out, Vector *max);
-  void CenterOrigin(ofstream *out, Vector *max);
-  void CenterGravity(ofstream *out, Vector *max);
+  void CenterOrigin(ofstream *out);
+  void CenterPeriodic(ofstream *out);
+  void CenterAtVector(ofstream *out, Vector *newcenter);
   void Translate(const Vector *x);
   void TranslatePeriodically(const Vector *trans);
@@ -261 +270 @@
   void Align(Vector *n);
   void Scale(double **factor);
-  void DetermineCenter(Vector &center);
+  void DeterminePeriodicCenter(Vector &center);
   Vector * DetermineCenterOfGravity(ofstream *out);
   Vector * DetermineCenterOfAll(ofstream *out);
-  void SetNameFromFilename(char *filename);
+  void SetNameFromFilename(const char *filename);
   void SetBoxDimension(Vector *dim);
   double * ReturnFullMatrixforSymmetric(double *cell_size);
   void ScanForPeriodicCorrection(ofstream *out);
+  bool VerletForceIntegration(ofstream *out, char *file, config &configuration);
+  void Thermostats(config &configuration, double ActualTemp, int Thermostat);
   void PrincipalAxisSystem(ofstream *out, bool DoRotate);
   double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
   Vector* FindEmbeddingHole(ofstream *out, molecule *srcmol);
 
-  bool VerletForceIntegration(char *file, double delta_t, bool IsAngstroem);
-
+
+  double ConstrainedPotential(ofstream *out, atom **permutation, int start, int end, double *constants, bool IsAngstroem);
+  double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
+  void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
+  bool LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration); 
+        
   bool CheckBounds(const Vector *x) const;
   void GetAlignvector(struct lsq_params * par) const;
@@ -349 +364 @@
   bool AddHydrogenCorrection(ofstream *out, char *path);
   bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
-  bool insert(molecule *mol);
+  void insert(molecule *mol);
   molecule * ReturnIndex(int index);
   bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
@@ -429 +444 @@
     char *databasepath;
 
+    int DoConstrainedMD;
+    int MaxOuterStep;
+    int Thermostat;
+    int *ThermostatImplemented;
+    char **ThermostatNames;
+    double TempFrequency;
+    double alpha;
+    double HooverMass;
+    double TargetTemp;
+    int ScaleTempStep;
+    
   private:
     char *mainname;
@@ -449 +475 @@
     int Seed;
 
-    int MaxOuterStep;
     int OutVisStep;
     int OutSrcStep;
-    double TargetTemp;
-    int ScaleTempStep;
     int MaxPsiStep;
     double EpsWannier;
@@ -501 +524 @@
   char *GetDefaultPath() const;
   void SetDefaultPath(const char *path);
+  void InitThermostats(ifstream *source);
 };
 

src/parser.cpp

@@ -56 +56 @@
 MatrixContainer::MatrixContainer() {
   Indices = NULL;
-  Header = (char *) Malloc(sizeof(char)*1023, "MatrixContainer::MatrixContainer: *Header");
+  Header = (char **) Malloc(sizeof(char)*1, "MatrixContainer::MatrixContainer: **Header");
   Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
   RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter");
+  ColumnCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *ColumnCounter");
+  Header[0] = NULL;
   Matrix[0] = NULL;
   RowCounter[0] = -1;
   MatrixCounter = 0;
-  ColumnCounter = 0;
+  ColumnCounter[0] = -1;
 };
 
@@ -70 +72 @@
   if (Matrix != NULL) {
     for(int i=MatrixCounter;i--;) {
-      if (RowCounter != NULL) {
+      if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
           for(int j=RowCounter[i]+1;j--;)
             Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
@@ -76 +78 @@
       }
     }
-    if ((RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
+    if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
       for(int j=RowCounter[MatrixCounter]+1;j--;)
         Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
@@ -87 +89 @@
       Free((void **)&Indices[i], "MatrixContainer::~MatrixContainer: *Indices[]");
     }
-  Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices");
-
-  Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
+  Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices"); 
+  
+  if (Header != NULL)
+    for(int i=MatrixCounter+1;i--;)
+      Free((void **)&Header[i], "MatrixContainer::~MatrixContainer: *Header[]");
+  Free((void **)&Header, "MatrixContainer::~MatrixContainer: **Header");
   Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
-};
-
+  Free((void **)&ColumnCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
+};
+
+/** Either copies index matrix from another MatrixContainer or initialises with trivial mapping if NULL.
+ * This either copies the index matrix or just maps 1 to 1, 2 to 2 and so on for all fragments.
+ * \param *Matrix pointer to other MatrixContainer
+ * \return true - copy/initialisation sucessful, false - dimension false for copying
+ */
+bool MatrixContainer::InitialiseIndices(class MatrixContainer *Matrix) 
+{
+  cout << "Initialising indices";
+  if (Matrix == NULL) {
+    cout << " with trivial mapping." << endl;
+    Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "MatrixContainer::InitialiseIndices: **Indices");
+    for(int i=MatrixCounter+1;i--;) {
+      Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
+      for(int j=RowCounter[i];j--;)
+        Indices[i][j] = j;
+    }
+  } else {
+    cout << " from other MatrixContainer." << endl;
+    if (MatrixCounter != Matrix->MatrixCounter)
+      return false;
+    Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "MatrixContainer::InitialiseIndices: **Indices");
+    for(int i=MatrixCounter+1;i--;) {
+      if (RowCounter[i] != Matrix->RowCounter[i])
+        return false;
+      Indices[i] = (int *) Malloc(sizeof(int)*Matrix->RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
+      for(int j=Matrix->RowCounter[i];j--;) {
+        Indices[i][j] = Matrix->Indices[i][j];
+        //cout << Indices[i][j] << "\t";
+      }
+      //cout << endl;
+    }
+  }
+  return true; 
+};
 
 /** Parsing a number of matrices.
@@ -121 +161 @@
   }
 
-  // skip some initial lines
+  // parse header
+  Header[MatrixNr] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseMatrix: *Header[]");
   for (int m=skiplines+1;m--;)
-    input.getline(Header, 1023);
-
+    input.getline(Header[MatrixNr], 1023);
+  
   // scan header for number of columns
-  line.str(Header);
+  line.str(Header[MatrixNr]);
   for(int k=skipcolumns;k--;)
-    line >> Header;
+    line >> Header[MatrixNr];
   //cout << line.str() << endl;
-  ColumnCounter=0;
+  ColumnCounter[MatrixNr]=0;
   while ( getline(line,token, '\t') ) {
     if (token.length() > 0)
-      ColumnCounter++;
+      ColumnCounter[MatrixNr]++;
   }
   //cout << line.str() << endl;
-  //cout << "ColumnCounter: " << ColumnCounter << "." << endl;
-  if (ColumnCounter == 0)
-    cerr << "ColumnCounter: " << ColumnCounter << " from file " << name << ", this is probably an error!" << endl;
-
+  //cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
+  if (ColumnCounter[MatrixNr] == 0)
+    cerr << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
+  
   // scan rest for number of rows/lines
   RowCounter[MatrixNr]=-1;    // counts one line too much
@@ -155 +196 @@
 
   // allocate matrix if it's not zero dimension in one direction
-  if ((ColumnCounter > 0) && (RowCounter[MatrixNr] > -1)) {
-    Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
-
+  if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
+    Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseMatrix: **Matrix[]");
+  
     // parse in each entry for this matrix
     input.clear();
     input.seekg(ios::beg);
     for (int m=skiplines+1;m--;)
-      input.getline(Header, 1023);    // skip header
-    line.str(Header);
+      input.getline(Header[MatrixNr], 1023);    // skip header
+    line.str(Header[MatrixNr]);
     for(int k=skipcolumns;k--;)  // skip columns in header too
       line >> filename;
-    strncpy(Header, line.str().c_str(), 1023);
+    strncpy(Header[MatrixNr], line.str().c_str(), 1023);  
     for(int j=0;j<RowCounter[MatrixNr];j++) {
-      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]");
       input.getline(filename, 1023);
       stringstream lines(filename);
@@ -174 +215 @@
       for(int k=skipcolumns;k--;)
         lines >> filename;
-      for(int k=0;(k<ColumnCounter) && (!lines.eof());k++) {
+      for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
         lines >> Matrix[MatrixNr][j][k];
         //cout << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
       }
       //cout << endl;
-      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
-      for(int j=ColumnCounter;j--;)
+      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]");
+      for(int j=ColumnCounter[MatrixNr];j--;)
         Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
     }
   } else {
-    cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl;
+    cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl; 
   }
   input.close();
@@ -233 +274 @@
 
   cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
+  Header = (char **) ReAlloc(Header, sizeof(char *)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: **Header"); // one more each for the total molecule
   Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
   RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
+  ColumnCounter = (int *) ReAlloc(ColumnCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *ColumnCounter");
   for(int i=MatrixCounter+1;i--;) {
     Matrix[i] = NULL;
+    Header[i] = NULL;
     RowCounter[i] = -1;
+    ColumnCounter[i] = -1;
   }
   for(int i=0; i < MatrixCounter;i++) {
@@ -252 +297 @@
 
 /** Allocates and resets the memory for a number \a MCounter of matrices.
- * \param *GivenHeader Header line
+ * \param **GivenHeader Header line for each matrix
  * \param MCounter number of matrices
  * \param *RCounter number of rows for each matrix
- * \param CCounter number of columns for all matrices
- * \return Allocation successful
- */
-bool MatrixContainer::AllocateMatrix(const char *GivenHeader, int MCounter, int *RCounter, int CCounter)
-{
-  Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader");
-  strncpy(Header, GivenHeader, 1023);
-
+ * \param *CCounter number of columns for each matrix
+ * \return Allocation successful 
+ */
+bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
+{
   MatrixCounter = MCounter;
-  ColumnCounter = CCounter;
-  Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
-  RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
+  Header = (char **) Malloc(sizeof(char *)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *Header");
+  Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: ***Matrix"); // one more each for the total molecule
+  RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *RowCounter");
+  ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *ColumnCounter");
   for(int i=MatrixCounter+1;i--;) {
+    Header[i] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::AllocateMatrix: *Header[i]");
+    strncpy(Header[i], GivenHeader[i], 1023);
     RowCounter[i] = RCounter[i];
-    Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
+    ColumnCounter[i] = CCounter[i];
+    Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::AllocateMatrix: **Matrix[]");  
     for(int j=RowCounter[i]+1;j--;) {
-      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
-      for(int k=ColumnCounter;k--;)
+      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]");
+      for(int k=ColumnCounter[i];k--;)
         Matrix[i][j][k] = 0.;
     }
@@ -286 +332 @@
   for(int i=MatrixCounter+1;i--;)
     for(int j=RowCounter[i]+1;j--;)
-      for(int k=ColumnCounter;k--;)
+      for(int k=ColumnCounter[i];k--;)
         Matrix[i][j][k] = 0.;
    return true;
@@ -299 +345 @@
   for(int i=MatrixCounter+1;i--;)
     for(int j=RowCounter[i]+1;j--;)
-      for(int k=ColumnCounter;k--;)
+      for(int k=ColumnCounter[i];k--;)
         if (fabs(Matrix[i][j][k]) > max)
           max = fabs(Matrix[i][j][k]);
@@ -315 +361 @@
   for(int i=MatrixCounter+1;i--;)
     for(int j=RowCounter[i]+1;j--;)
-      for(int k=ColumnCounter;k--;)
+      for(int k=ColumnCounter[i];k--;)
         if (fabs(Matrix[i][j][k]) < min)
           min = fabs(Matrix[i][j][k]);
@@ -331 +377 @@
 {
   for(int j=RowCounter[MatrixCounter]+1;j--;)
-    for(int k=skipcolumns;k<ColumnCounter;k++)
+    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
       Matrix[MatrixCounter][j][k] = value;
    return true;
@@ -344 +390 @@
 {
   for(int j=RowCounter[MatrixCounter]+1;j--;)
-    for(int k=skipcolumns;k<ColumnCounter;k++)
+    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
       Matrix[MatrixCounter][j][k] = values[j][k];
    return true;
 };
 
-/** Sums the energy with each factor and put into last element of \a ***Energies.
+/** Sums the entries with each factor and put into last element of \a ***Matrix.
  * Sums over "E"-terms to create the "F"-terms
- * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
+ * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
  * \param Order bond order
@@ -384 +430 @@
               }
               if (Order == SubOrder) { // equal order is always copy from Energies
-                for(int l=ColumnCounter;l--;) // then adds/subtract each column
+                for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column 
                   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
               } else {
-                for(int l=ColumnCounter;l--;)
+                for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;)
                   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
               }
             }
-            //if ((ColumnCounter>1) && (RowCounter[0]-1 >= 1))
+            //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
              //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
           }
@@ -426 +472 @@
       return false;
     }
-    output << Header << endl;
+    output << Header[i] << endl;
     for(int j=0;j<RowCounter[i];j++) {
-      for(int k=0;k<ColumnCounter;k++)
+      for(int k=0;k<ColumnCounter[i];k++)
         output << scientific << Matrix[i][j][k] << "\t";
       output << endl;
@@ -455 +501 @@
     return false;
   }
-  output << Header << endl;
+  output << Header[MatrixCounter] << endl;
   for(int j=0;j<RowCounter[MatrixCounter];j++) {
-    for(int k=0;k<ColumnCounter;k++)
+    for(int k=0;k<ColumnCounter[MatrixCounter];k++)
       output << scientific << Matrix[MatrixCounter][j][k] << "\t";
     output << endl;
@@ -485 +531 @@
 /** Sums the energy with each factor and put into last element of \a EnergyMatrix::Matrix.
  * Sums over the "F"-terms in ANOVA decomposition.
- * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
+ * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
  * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
@@ -498 +544 @@
     for(int i=KeySet.FragmentsPerOrder[Order];i--;)
       for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
-        for(int k=ColumnCounter;k--;)
+        for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
           Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
   else
     for(int i=KeySet.FragmentsPerOrder[Order];i--;)
       for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
-        for(int k=ColumnCounter;k--;)
+        for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
           Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
   return true;
@@ -524 +570 @@
     // count maximum of columns
     RowCounter[MatrixCounter] = 0;
-    for(int j=0; j < MatrixCounter;j++) // (energy matrix might be bigger than number of atoms in terms of rows)
+    ColumnCounter[MatrixCounter] = 0;
+    for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
       if (RowCounter[j] > RowCounter[MatrixCounter])
         RowCounter[MatrixCounter] = RowCounter[j];
+      if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
+        ColumnCounter[MatrixCounter] = ColumnCounter[j];
+    }
     // allocate last plus one matrix
-    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
+    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
     Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     for(int j=0;j<=RowCounter[MatrixCounter];j++)
-      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
-
+      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+    
     // try independently to parse global energysuffix file if present
     strncpy(filename, name, 1023);
@@ -589 +639 @@
 
 /** Sums the forces and puts into last element of \a ForceMatrix::Matrix.
- * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
+ * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
  * \param Order bond order
@@ -609 +659 @@
       if (j != -1) {
         //if (j == 0) cout << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
-        for(int k=2;k<ColumnCounter;k++)
+        for(int k=2;k<ColumnCounter[MatrixCounter];k++)
           Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
       }
@@ -656 +706 @@
     input.close();
 
+    ColumnCounter[MatrixCounter] = 0;
+    for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
+      if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
+        ColumnCounter[MatrixCounter] = ColumnCounter[j];
+    }
+  
     // allocate last plus one matrix
-    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
+    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
     Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     for(int j=0;j<=RowCounter[MatrixCounter];j++)
-      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
+
+    // try independently to parse global forcesuffix file if present
+    strncpy(filename, name, 1023);
+    strncat(filename, prefix, 1023-strlen(filename));
+    strncat(filename, suffix.c_str(), 1023-strlen(filename));
+    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
+  }
+ 
+
+  return status;
+};
+
+// ======================================= CLASS HessianMatrix =============================
+
+/** Parsing force Indices of each fragment
+ * \param *name directory with \a ForcesFile
+ * \return parsing successful 
+ */
+bool HessianMatrix::ParseIndices(char *name) 
+{
+  ifstream input;
+  char *FragmentNumber = NULL;
+  char filename[1023];
+  stringstream line;
+  
+  cout << "Parsing hessian indices for " << MatrixCounter << " matrices." << endl;
+  Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "HessianMatrix::ParseIndices: **Indices");
+  line << name << FRAGMENTPREFIX << FORCESFILE;
+  input.open(line.str().c_str(), ios::in);
+  //cout << "Opening " << line.str() << " ... "  << input << endl;
+  if (input == NULL) {
+    cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
+    return false;
+  }
+  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
+    // get the number of atoms for this fragment
+    input.getline(filename, 1023);
+    line.str(filename);
+    // parse the values
+    Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "HessianMatrix::ParseIndices: *Indices[]");
+    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
+    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
+    Free((void **)&FragmentNumber, "HessianMatrix::ParseIndices: *FragmentNumber");
+    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
+      line >> Indices[i][j];
+      //cout << " " << Indices[i][j];
+    }
+    //cout << endl;
+  }
+  Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "HessianMatrix::ParseIndices: *Indices[]");
+  for(int j=RowCounter[MatrixCounter];j--;) {
+    Indices[MatrixCounter][j] = j;
+  }
+  input.close();
+  return true;
+};
+
+
+/** Sums the hessian entries and puts into last element of \a HessianMatrix::Matrix.
+ * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
+ * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
+ * \param Order bond order
+ *  \param sign +1 or -1
+ * \return true if summing was successful
+ */
+bool HessianMatrix::SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
+{
+  int FragmentNr;
+  // sum forces
+  for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
+    FragmentNr = KeySet.OrderSet[Order][i];
+    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
+      int j = Indices[ FragmentNr ][l];
+      if (j > RowCounter[MatrixCounter]) {
+        cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
+        return false;
+      }
+      if (j != -1) {
+        for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
+          int k = Indices[ FragmentNr ][m];
+          if (k > ColumnCounter[MatrixCounter]) {
+            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
+            return false;
+          }
+          if (k != -1) {
+            //cout << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
+            Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
+          }
+        }
+      }
+    }
+  }
+  return true;
+};
+
+/** Constructor for class HessianMatrix.
+ */
+HessianMatrix::HessianMatrix() : MatrixContainer() 
+{
+   IsSymmetric = true;
+}
+
+/** Sums the hessian entries with each factor and put into last element of \a ***Matrix.
+ * Sums over "E"-terms to create the "F"-terms
+ * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
+ * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
+ * \param Order bond order
+ * \return true if summing was successful
+ */
+bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
+{
+  // go through each order
+  for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
+    //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
+    // then go per order through each suborder and pick together all the terms that contain this fragment
+    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
+      for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
+        if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
+          //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
+          // if the fragment's indices are all in the current fragment
+          for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
+            int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
+            //cout << "Current row index is " << k << "/" << m << "." << endl;
+            if (m != -1) { // if it's not an added hydrogen
+              for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
+                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
+                if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
+                  m = l;
+                  break;  
+                }
+              }
+              //cout << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
+              if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
+                cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
+                return false;
+              }
+              
+              for(int l=0;l<ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l++) { 
+                int n = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][l];
+                //cout << "Current column index is " << l << "/" << n << "." << endl;
+                if (n != -1) { // if it's not an added hydrogen
+                  for (int p=0;p<ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
+                    //cout << "Comparing " << n << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
+                    if (n == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p]) {
+                      n = p;
+                      break;  
+                    }
+                  }
+                  //cout << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
+                  if (n > ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
+                    cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
+                    return false;
+                  }
+                  if (Order == SubOrder) { // equal order is always copy from Energies
+                    //cout << "Adding " << MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
+                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
+                  } else {
+                    //cout << "Subtracting " << Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
+                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
+                  }
+                }
+              }
+            }
+            //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
+             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
+          }
+        } else {
+          //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
+        }
+      }
+    }
+   //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
+  }
+  
+  return true;
+};
+
+/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
+ * \param *name directory with files
+ * \param *prefix prefix of each matrix file
+ * \param *suffix suffix of each matrix file
+ * \param skiplines number of inital lines to skip
+ * \param skiplines number of inital columns to skip
+ * \return parsing successful
+ */ 
+bool HessianMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
+{
+  char filename[1023];
+  ifstream input;
+  stringstream file;
+  int nr;
+  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
+
+  if (status) {
+    // count number of atoms for last plus one matrix
+    file << name << FRAGMENTPREFIX << KEYSETFILE;
+    input.open(file.str().c_str(), ios::in);
+    if (input == NULL) {
+      cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
+      return false;
+    }
+    RowCounter[MatrixCounter] = 0;
+    ColumnCounter[MatrixCounter] = 0;
+    while (!input.eof()) {
+      input.getline(filename, 1023);
+      stringstream zeile(filename);
+      while (!zeile.eof()) {
+        zeile >> nr;
+        //cout << "Current index: " << nr << "." << endl;
+        if (nr > RowCounter[MatrixCounter]) {
+          RowCounter[MatrixCounter] = nr;
+          ColumnCounter[MatrixCounter] = nr;
+        }
+      }
+    }
+    RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
+    ColumnCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
+    input.close();
+  
+    // allocate last plus one matrix
+    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
+    Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
+    for(int j=0;j<=RowCounter[MatrixCounter];j++)
+      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
 
     // try independently to parse global forcesuffix file if present

src/parser.hpp

@@ -19 +19 @@
 
 #define EnergySuffix ".energy.all"
+#define EnergyFragmentSuffix ".energyfragment.all"
+#define ForcesSuffix ".forces.all"
+#define ForceFragmentSuffix ".forcefragment.all"
 #define HcorrectionSuffix ".Hcorrection.all"
-#define ForcesSuffix ".forces.all"
+#define HcorrectionFragmentSuffix ".Hcorrectionfragment.all"
+#define HessianSuffix ".hessian_xx.all"
+#define HessianFragmentSuffix ".hessianfragment_xx.all"
+#define OrderSuffix ".Order"
 #define ShieldingSuffix ".sigma_all.csv"
 #define ShieldingPASSuffix ".sigma_all_PAS.csv"
 #define ShieldingFragmentSuffix ".sigma_all_fragment.all"
 #define ShieldingPASFragmentSuffix ".sigma_all_PAS_fragment.all"
+#define ChiSuffix ".chi_all.csv"
+#define ChiPASSuffix ".chi_all_PAS.csv"
+#define ChiFragmentSuffix ".chi_all_fragment.all"
+#define ChiPASFragmentSuffix ".chi_all_PAS_fragment.all"
 #define TimeSuffix ".speed"
-#define EnergyFragmentSuffix ".energyfragment.all"
-#define HcorrectionFragmentSuffix ".Hcorrectionfragment.all"
-#define ForceFragmentSuffix ".forcefragment.all"
-#define OrderSuffix ".Order"
 
 // ======================================= FUNCTIONS ==========================================
@@ -43 +49 @@
     double ***Matrix;
     int **Indices;
-    char *Header;
+    char **Header;
     int MatrixCounter;
     int *RowCounter;
-    int ColumnCounter;
-
+    int *ColumnCounter;
+  
   MatrixContainer();
-  virtual ~MatrixContainer();
-
+  ~MatrixContainer();
+  
+  bool InitialiseIndices(class MatrixContainer *Matrix = NULL);
   bool ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr);
   virtual bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns);
-  bool AllocateMatrix(const char *GivenHeader, int MCounter, int *RCounter, int CCounter);
+  bool AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter);
   bool ResetMatrix();
   double FindMinValue();
@@ -84 +91 @@
 };
 
+// ======================================= CLASS HessianMatrix =============================
+
+class HessianMatrix : public MatrixContainer { 
+  public: 
+    HessianMatrix();
+    //~HessianMatrix();
+    bool ParseIndices(char *name);
+    bool SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order);
+    bool SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign);
+    bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns);
+  private:
+    bool IsSymmetric;
+};
+
 // ======================================= CLASS KeySetsContainer =============================