Changeset 390248

Timestamp:

Jul 24, 2008, 2:00:19 PM (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:

958457

Parents:

c685eb

Message:

BIG change: Hcorrection included and bugfix of analyzer (wrt forces)

• Hcorrection is supplied via molecuilder which gives correction values if hydrogen atoms are to close and thus already feel each other's influence. This correction energy matrix is also parsed ans subtracted from the approximated energy matrix generated from the fragments. This impacts both on joiner, analyzer and molecuilder
• Forces were not working which had multiple causes:
  • first, we stepped back down to absolute errors which rather gives a feel about convergence (1e-6 is simply small and neglegible)
  • there may have been bugs where (either in Force or in ForceFragments) adding up took place, this is cleaned up
  • more routines have been abstracted from analyzer into datacreator functions
    • the order of ions had changed, re-run with current molecuilder, calculation and subsequent joining/analyzing brought correct results finally
• plots were bad still due to wrong axes (two more functions finding max/min values), use of wrong columns in plot file

Location:

src

Files:

• analyzer.cpp (modified) (12 diffs)
• datacreator.cpp (modified) (9 diffs)
• datacreator.hpp (modified) (2 diffs)
• defs.hpp (modified) (1 diff)
• joiner.cpp (modified) (8 diffs)
• moleculelist.cpp (modified) (3 diffs)
• molecules.cpp (modified) (1 diff)
• molecules.hpp (modified) (1 diff)
• parser.cpp (modified) (8 diffs)
• parser.hpp (modified) (5 diffs)

src/analyzer.cpp

@@ -25 +25 @@
   periodentafel *periode = NULL; // and a period table of all elements
   EnergyMatrix Energy;
+  EnergyMatrix Hcorrection;
   ForceMatrix Force;
   ForceMatrix Shielding;
@@ -30 +31 @@
   MatrixContainer Time;
   EnergyMatrix EnergyFragments;
+  EnergyMatrix HcorrectionFragments;
   ForceMatrix ForceFragments;
   ForceMatrix ShieldingFragments;
@@ -45 +47 @@
   stringstream Orderxrange;
   stringstream Fragmentxrange;
-
+  stringstream yrange;
+  char *dir = NULL;
+  bool Hcorrected = true;
+  double norm;
+  
   cout << "ANOVA Analyzer" << endl;
   cout << "==============" << endl;
@@ -57 +63 @@
     cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl;
     return 1;
-  }
+  } else {
+    dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir");
+    strcpy(dir, "/");
+    strcat(dir, argv[2]);
+  }
+  
   if (argc > 4) {
     cout << "Loading periodentafel." << endl;
@@ -71 +82 @@
   
   // ------------- Parse through all Fragment subdirs --------
-  if (!Energy.ParseMatrix(argv[1], argv[2], EnergySuffix,0,0)) return 1;
-  if (!Force.ParseMatrix(argv[1], argv[2], ForcesSuffix,0,0)) return 1;
-  if (!Time.ParseMatrix(argv[1], argv[2], TimeSuffix, 10,1)) return 1;
+  if (!Energy.ParseMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
+  Hcorrected = Hcorrection.ParseMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
+  if (!Force.ParseMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
+  //if (!Time.ParseMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if (!Shielding.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
-    if (!ShieldingPAS.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
+    if (!Shielding.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
+    if (!ShieldingPAS.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
   }
 
   // ---------- Parse the TE Factors into an array -----------------
   if (!Energy.ParseIndices()) return 1;
+  if (Hcorrected) 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;
 
   // ---------- Parse the shielding indices into an array ---------------
@@ -96 +111 @@
   
   // ---------- Parse fragment files created by 'joiner' into an array -------------
-  if (!EnergyFragments.ParseMatrix(argv[1], argv[2], EnergyFragmentSuffix,0,0)) return 1;
-  if (!ForceFragments.ParseMatrix(argv[1], argv[2], ForceFragmentSuffix,0,0)) return 1;
+  if (!EnergyFragments.ParseMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
+  if (Hcorrected) HcorrectionFragments.ParseMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
+  if (!ForceFragments.ParseMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if (!ShieldingFragments.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
-    if (!ShieldingPASFragments.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
+    if (!ShieldingFragments.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
+    if (!ShieldingPASFragments.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
   }
 
@@ -143 +159 @@
   }
   
-  output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
-  Time.SetLastMatrix(0., 0);
-  for (int 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];
-        }
-  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 Times" << endl << "===============" << endl << Time.Header << endl;
+//  Time.SetLastMatrix(0., 0);
+//  for (int 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];
+//        }
+//  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.close();
 
@@ -165 +181 @@
   // ======================================= Creating the data files ==============================================================
 
-  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
-  if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
-  Time.SetLastMatrix(0., 0);
-  output << "#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];
-        }
-    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
-    for(int k=0;k<Time.ColumnCounter;k++)
-      output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][1][k];
-    output << endl;
-  }
-  output.close();
+//  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
+//  if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
+//  Time.SetLastMatrix(0., 0);
+//  output << "#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];
+//        }
+//    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+//    for(int k=0;k<Time.ColumnCounter;k++)
+//      output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][1][k];
+//    output << endl;
+//  }
+//  output.close();
 
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
-  if (!EnergyFragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0)) return 1;
-  if (!CreateDataEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "DeltaEnergies-Order", "Plot of error between approximated and full energies energies versus the Bond Order", datum)) return 1;
+  if (!CreateDataDeltaEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "DeltaEnergies-Order", "Plot of error between approximated and full energies energies versus the Bond Order", datum)) return 1;
   
   // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
-  if (!EnergyFragments.SetLastMatrix(1.,0)) return 1;
-  if (!CreateDataEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
+  if (!CreateDataEnergyOrder(EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
   
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
-  if (!ForceFragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0)) return 1;
-  if (!OpenOutputFile(output, argv[3], "DeltaForces-Order.dat" )) return false;
-  if (!OpenOutputFile(output2, argv[3], "DeltaMinForces-Order.dat" )) return false;
-  if (!OpenOutputFile(output3, argv[3], "DeltaMeanForces-Order.dat" )) return false;
-  if (!OpenOutputFile(output4, argv[3], "DeltaMaxForces-Order.dat" )) return false;
-  cout << "Plot of per atom and min/mean/max error between approximated forces and full forces versus the Bond Order ... " << endl;
-  output << "# Plot of error between approximated forces and full forces versus the Bond Order, created on " << datum;
-  output << "# AtomNo" << Force.Header << endl;
-  output2 << "# Plot of min error between approximated forces and full forces versus the Bond Order, created on " << datum;
-  output2 << "# Order\tFrag.No.\t" << Force.Header << endl;
-  output3 << "# Plot of mean error between approximated forces and full forces versus the Bond Order, created on " << datum;
-  output3 << "# Order\tFrag.No.\t" << Force.Header << endl;
-  output4 << "# Plot of max error between approximated forces and full forces versus the Bond Order, created on " << datum;
-  output4 << "# Order\tFrag.No.\t" << Force.Header << endl;
-  int FragmentCounter = 0;
-  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
-    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
-      for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
-        int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
-        if (j > Force.RowCounter[Force.MatrixCounter]) {
-          cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
-          return 1;
-        }
-        if (j != -1)
-          for(int k=0;k<Force.ColumnCounter;k++) {
-            Force.Matrix[Force.MatrixCounter][j][k] -= ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
-          }
-      }
-      FragmentCounter++;
-    }
-    // errors per atom
-    output << "#Order\t" << BondOrder+1 << endl;
-    for(int j=0;j<Force.RowCounter[ Force.MatrixCounter ];j++) { 
-      output << Force.Indices[Force.MatrixCounter][j] << "\t";
-      for (int l=0;l<Force.ColumnCounter;l++)
-        if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ j ][l]) < MYEPSILON)
-          output << scientific << Force.Matrix[Force.MatrixCounter][ j ][l] << "\t";
-        else
-          output << scientific << (Force.Matrix[Force.MatrixCounter][ j ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ j ][l]) << "\t";
-      output << endl;
-    }
-    output << endl;
-    
-    // min error
-    output2 << BondOrder+1 << "\t" << FragmentCounter;
-    CreateMinimumForce(Force, Force.MatrixCounter);
-    for (int l=0;l<Force.ColumnCounter;l++)
-      if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
-        output2 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
-      else
-        output2 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
-    output2 << endl;
-    
-    // mean error
-    output3 << BondOrder+1 << "\t" << FragmentCounter;
-    CreateMeanForce(Force, Force.MatrixCounter);
-    for (int l=0;l<Force.ColumnCounter;l++)
-      if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
-        output3 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
-      else
-        output3 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
-    output3 << endl;
-    
-    // max error
-    output4 << BondOrder+1 << "\t" << FragmentCounter;
-    CreateMaximumForce(Force, Force.MatrixCounter);
-    for (int l=0;l<Force.ColumnCounter;l++)
-      if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
-        output4 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
-      else
-        output4 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
-    output4 << endl;
-  }
-  output.close();
-  output2.close();
-  output3.close();
-  output4.close();
+  if (!CreateDataDeltaForcesOrderPerAtom(Force, ForceFragments, KeySet, argv[3], "DeltaForces-Order", "Plot of error between approximated forces and full forces versus the Bond Order", datum)) return 1;
+
+  // min force
+  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMinForces-Order", "Plot of min error between approximated forces and full forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
+
+  // mean force
+  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMeanForces-Order", "Plot of mean error between approximated forces and full forces versus the Bond Order", datum, CreateMeanForce)) return 1;
+
+  // max force
+  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMaxForces-Order", "Plot of max error between approximated forces and full forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
 
   // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
-  if(!OpenOutputFile(output, argv[3], "Forces-Order.dat")) return 1;
-  cout << "Plot of approximated forces versus the Bond Order ... " << endl;
-  output << "# Plot of approximated forces versus the Bond Order, created on " << datum;
-  output << "# AtomNo" << Force.Header << endl;
-  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
-    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
-      for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
-        int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
-        if (j > Force.RowCounter[Force.MatrixCounter]) {
-          cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
-          return 1;
-        }
-        if (j != -1)
-          for(int k=0;k<Force.ColumnCounter;k++) {
-            Force.Matrix[Force.MatrixCounter][j][k] -= ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
-          }
-      }
-    }
-    // errors per atom
-    output << "#Order\t" << BondOrder+1 << endl;
-    for(int j=0;j<Force.RowCounter[ Force.MatrixCounter ];j++) { 
-      output << Force.Indices[Force.MatrixCounter][j] << "\t";
-      for (int l=0;l<Force.ColumnCounter;l++)
-         output << Force.Matrix[Force.MatrixCounter][ j ][l] << "\t";
-      output << endl;
-    }
-    output << endl;
-  }
-  output.close();
-
-  // min force
-  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
-
-  // mean force
-  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;
-
-  // max force
-  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
+  if (!CreateDataForcesOrderPerAtom(ForceFragments, KeySet, argv[3], "Forces-Order", "Plot of approximated forces versus the Bond Order", datum)) return 1;
+
+  // min force
+  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
+
+  // mean force
+  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;
+
+  // max force
+  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
 
   // ++++++++++++++++++++++++++++++++++++++Plotting vector sum (should be 0) vs. bond order
-  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;
+  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;
 
   // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
@@ -342 +261 @@
   Orderxrange << "[1:" << KeySet.Order << "]";
   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 (!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(), "[1e-8:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
+  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;
   
   // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
@@ -353 +273 @@
 
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
-  // min force
-  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMinForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in min force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
-
-  // mean force
-  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMeanForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in mean force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;
-
-  // max force
-  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMaxForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in max force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
+  yrange.str("[1e-8:1e+0]");
+  // min force
+  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMinForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in min force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
+
+  // mean force
+  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMeanForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in mean force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;
+
+  // max force
+  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMaxForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in max force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
 
   // min/mean/max comparison for total force
@@ -395 +316 @@
 
   // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
-  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-Fragment", 5, "below", "y", "", 1, 5, "fragment number", "Energies of each fragment [Ht]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with points", AbsEnergyPlotLine)) return 1;
-  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Energies of each fragment [Ht]", Orderxrange.str().c_str(), "[1e-10:1e+1]", "1" , "with points", AbsEnergyPlotLine)) return 1;
-  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Maximum fragment energy [Ht]", Orderxrange.str().c_str(), "[1e-6:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
-  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Minimum fragment energy [Ht]", Orderxrange.str().c_str(), "[1e-6:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
+  yrange.str("");
+  yrange << "[" << EnergyFragments.FindMinValue() << ":" << EnergyFragments.FindMaxValue() << "]";
+  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-Fragment", 5, "below", "y", "", 1, 5, "fragment number", "Energies of each fragment [Ht]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with points", AbsEnergyPlotLine)) return 1;
+  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Energies of each fragment [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with points", AbsEnergyPlotLine)) return 1;
+  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Maximum fragment energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
+  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Minimum fragment energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
 
   // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment
+  yrange.str("");
+  yrange << "[" << ForceFragments.FindMinValue() << ":" << ForceFragments.FindMaxValue()<< "]";
   // min
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "minimum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "minimum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
     
   // mean
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "mean of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "mean of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
     
   // max
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "maximum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "maximum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
 
   // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment per bond order
   // min
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "minimum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "minimum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
     
   // mean
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "mean of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "mean of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
     
   // max
-  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "maximum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
+  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "maximum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
  
   // create Makefile
@@ -434 +359 @@
   // ++++++++++++++++ exit ++++++++++++++++++++++++++++++++++
   delete(periode);
+  Free((void **)&dir, "main: *dir");
   cout << "done." << endl;
   return 0;

src/datacreator.cpp

@@ -30 +30 @@
 
 /** Plots an energy vs. order.
- * \param Energy EnergyMatrix class containing matrix values (last matrix is created by summing over Fragments)
- * \param EnergyFragments EnergyMatrix class containing matrix values
+ * \param &Fragments EnergyMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
@@ -37 +36 @@
  * \return true if file was written successfully
  */
-bool CreateDataEnergyOrder(class MatrixContainer &Energy, class MatrixContainer &EnergyFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
+bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
 {
   stringstream filename;
@@ -46 +45 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "#Order\tFrag.No.\t" << Energy.Header << endl;
+  output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
-      for(int j=Energy.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
-        for(int k=Energy.ColumnCounter;k--;)
-          Energy.Matrix[Energy.MatrixCounter][j][k] -= EnergyFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+      for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
+        for(int k=Fragments.ColumnCounter;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<Energy.ColumnCounter;l++)
-      if (fabs(EnergyFragments.Matrix[EnergyFragments.MatrixCounter][ EnergyFragments.RowCounter[EnergyFragments.MatrixCounter]-1 ][l]) < MYEPSILON)
-        output << scientific << "\t" << Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l];
+    for (int l=0;l<Fragments.ColumnCounter;l++)
+      output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plots an energy error vs. order.
+ * \param &Energy EnergyMatrix class containing reference values (in MatrixCounter matrix)
+ * \param &Fragments EnergyMatrix 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
+ * \return true if file was written successfully
+ */
+bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, 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 << "#Order\tFrag.No.\t" << Fragments.Header << 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--;)
+          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++)
+      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];
       else
-        output << scientific << "\t" << (Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l] / EnergyFragments.Matrix[EnergyFragments.MatrixCounter][ EnergyFragments.RowCounter[EnergyFragments.MatrixCounter]-1 ][l]);
+        output << scientific << "\t" << (Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l] / Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]);
     output << endl;
   }
@@ -66 +99 @@
 
 /** Plot forces vs. order.
- */
-bool CreateDataForcesOrder(class MatrixContainer &Force, class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
+ * \param &Fragments ForceMatrix 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 CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
 {
   stringstream filename;
@@ -76 +115 @@
   cout << msg << endl;
   output << "# " << msg << ", created on " << datum;
-  output << "# Order\tFrag.No.\t" << Force.Header << endl;
-  for(int j=0;j<=Force.RowCounter[ Force.MatrixCounter ];j++)
-    for(int k=0;k<Force.ColumnCounter;k++)
-       Force.Matrix[Force.MatrixCounter][j][k] = 0.;
-  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
-    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
-      for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
-        int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
-        if (j > Force.RowCounter[Force.MatrixCounter]) {
-          cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
-          return 1;
-        }
-        if (j != -1)
-          for(int k=Force.ColumnCounter;k--;) {
-            Force.Matrix[Force.MatrixCounter][j][k] += ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
-          }
+  output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
+  Fragments.SetLastMatrix(0.,0);
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
+    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+    CreateForce(Fragments, Fragments.MatrixCounter);
+    for (int l=0;l<Fragments.ColumnCounter;l++)
+       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot forces error vs. order.
+ * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
+ * \param &Fragments ForceMatrix 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 CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
+{
+  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 << "# Order\tFrag.No.\t" << Fragments.Header << endl;
+  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
+    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+    CreateForce(Fragments, Fragments.MatrixCounter);
+    for (int l=0;l<Fragments.ColumnCounter;l++)
+       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot forces error vs. vs atom vs. order.
+ * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
+ * \param &Fragments ForceMatrix 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 CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
+{
+  stringstream filename;
+  ofstream output;
+  double norm = 0.;
+
+  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 << endl;
+  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    //cout << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
+    // errors per atom
+    output << "#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;l++) {
+        if (((l+1) % 3) == 0) {
+          norm = 0.;
+          for (int m=0;m<NDIM;m++)
+            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";
+//        else
+//          output << scientific << (Fragments.Matrix[Fragments.MatrixCounter][ j ][l] / norm) << "\t";
       }
-    }
-    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
-    CreateForce(Force, Force.MatrixCounter);
-    for (int l=0;l<Force.ColumnCounter;l++)
-       output << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
+      output << endl;
+    }
+    output << endl;
+  }
+  output.close();
+  return true;
+};
+
+/** Plot forces error vs. vs atom vs. order.
+ * \param &Fragments ForceMatrix 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 CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, 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 << endl;
+  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+    //cout << "Current order is " << BondOrder << "." << endl;
+    Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
+    // errors per atom
+    output << "#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;l++)
+        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
+      output << endl;
+    }
     output << endl;
   }
@@ -286 +426 @@
 };
 
+/** Leaves the Force.Matrix as it is.
+ * \param Force ForceMatrix class containing matrix values
+ * \param MatrixNumber the index for the ForceMatrix::matrix array
+*/
+void CreateSameForce(class MatrixContainer &Force, int MatrixNumber)
+{
+  // does nothing
+};
+
 /** Adds vectorwise all forces.
  * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
@@ -295 +444 @@
   for (int l=Force.ColumnCounter;l--;)
     Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
-  for (int l=5;l<Force.ColumnCounter;l++) {
+  for (int l=0;l<Force.ColumnCounter;l++) {
     for (int k=Force.RowCounter[MatrixNumber];k--;)
       Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] += Force.Matrix[MatrixNumber][k][l];
@@ -375 +524 @@
 
   line >> item;
-  for (int i=3; i< Energy.ColumnCounter;i++) {
-    line >> item;
-    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
-    if (i != (Energy.ColumnCounter-1))
+  for (int i=2; i<= Energy.ColumnCounter;i++) {
+    line >> item;
+    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(abs($" << i+2 << ")) " << uses;
+    if (i != (Energy.ColumnCounter))
       output << ", \\";
     output << endl;
@@ -397 +546 @@
 
   line >> item;
-  for (int i=3; i< Energy.ColumnCounter;i++) {
-    line >> item;
-    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":" << i+1 << " " << uses;
-    if (i != (Energy.ColumnCounter-1))
+  for (int i=2; i<= Energy.ColumnCounter;i++) {
+    line >> item;
+    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":" << i+2 << " " << uses;
+    if (i != (Energy.ColumnCounter))
       output << ", \\";
     output << endl;

src/datacreator.hpp

@@ -19 +19 @@
 bool OpenOutputFile(ofstream &output, const char *dir, const char *filename);
 
-bool CreateDataEnergyOrder(class MatrixContainer &Energy, class MatrixContainer &EnergyFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 
-bool CreateDataForcesOrder(class MatrixContainer &Force, class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int));
+bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 
+bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 
+bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int));
+bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum);
+bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int));
+bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 
 bool CreateDataFragment(class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int));
 bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int));
@@ -30 +34 @@
 void CreateMeanForce(class MatrixContainer &Force, int MatrixNumber);
 void CreateMaximumForce(class MatrixContainer &Force, int MatrixNumber);
+void CreateSameForce(class MatrixContainer &Force, int MatrixNumber);
 void CreateVectorSumForce(class MatrixContainer &Force, int MatrixNumber);
 

src/defs.hpp

@@ -49 +49 @@
 #define TEFACTORSFILE "TE-Factors.dat"    //!< default filename of BOSSANOVA total energy factors file
 #define FORCESFILE "Forces-Factors.dat"    //!< default filename of BOSSANOVA force factors file
+#define HCORRECTIONSUFFIX "Hcorrection.dat"    //!< default filename of BOSSANOVA H correction file (unwanted saturation interaction)
+#define FITCONSTANTSUFFIX "FitConstant.dat"   //!< suffix of default filename of BOSSANOVA fit constants file (unwanted saturation interaction)
 #define SHIELDINGSUFFIX "sigma_all.csv"                //!< default filename of BOSSANOVA shieldings file
 #define SHIELDINGPASSUFFIX "sigma_all_PAS.csv"                 //!< default filename of BOSSANOVA shieldings PAS file

src/joiner.cpp

@@ -19 +19 @@
   periodentafel *periode = NULL; // and a period table of all elements
   EnergyMatrix Energy;
+  EnergyMatrix Hcorrection;
   ForceMatrix Force;
   EnergyMatrix EnergyFragments;
+  EnergyMatrix HcorrectionFragments;
   ForceMatrix ForceFragments;
   ForceMatrix Shielding;
@@ -28 +30 @@
   KeySetsContainer KeySet;  
   stringstream prefix;
+  char *dir = NULL;
+  bool Hcorrected = true;
 
   cout << "Joiner" << endl;
@@ -39 +43 @@
     cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl;
     return 1;
+  } else {
+    dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir");
+    strcpy(dir, "/");
+    strcat(dir, argv[2]);
   }
   if (argc > 3) {
@@ -52 +60 @@
   
   // ------------- Parse through all Fragment subdirs --------
-  if (!Energy.ParseMatrix(argv[1], argv[2], EnergySuffix, 0,0)) return 1;
-  if (!Force.ParseMatrix(argv[1], argv[2], ForcesSuffix, 0,0)) return 1;
+  if (!Energy.ParseMatrix(argv[1], dir, EnergySuffix, 0,0)) return 1;
+  Hcorrected = Hcorrection.ParseMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0);
+  if (!Force.ParseMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if (!Shielding.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
-    if (!ShieldingPAS.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
+    if (!Shielding.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
+    if (!ShieldingPAS.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
   }
 
   // ---------- Parse the TE Factors into an array -----------------
   if (!Energy.ParseIndices()) return 1;
+  if (Hcorrected) Hcorrection.ParseIndices();
   
   // ---------- Parse the Force indices into an array ---------------
@@ -76 +86 @@
   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 (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
   if (periode != NULL) { // also look for PAS values
@@ -86 +97 @@
   if(!Force.SetLastMatrix(0., 2)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if(!Shielding.SetLastMatrix(0., 0)) return 1;
+    if(!Shielding.SetLastMatrix(0., 2)) return 1;
     if(!ShieldingPAS.SetLastMatrix(0., 2)) return 1;
   }
@@ -97 +108 @@
     cout << "Summing energy of order " << BondOrder+1 << " ..." << endl;
     if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return 1;
-    if (!Energy.SumSubEnergy(EnergyFragments, KeySet, BondOrder)) return 1;
+    if (Hcorrected) { 
+      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 
+      if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1;
     // --------- sum up Forces --------------------
     cout << "Summing forces of order " << BondOrder+1 << " ..." << endl;
     if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return 1;
-    if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder)) return 1;
+    if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return 1;
     if (periode != NULL) { // also look for PAS values
       cout << "Summing shieldings of order " << BondOrder+1 << " ..." << endl;
       if (!ShieldingFragments.SumSubManyBodyTerms(Shielding, KeySet, BondOrder)) return 1;
-      if (!Shielding.SumSubForces(ShieldingFragments, 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)) return 1;
+      if (!ShieldingPAS.SumSubForces(ShieldingPASFragments, KeySet, BondOrder, 1.)) return 1;
     }
 
     // --------- write the energy and forces file --------------------
     prefix.str(" ");
-    prefix << argv[2] << OrderSuffix << (BondOrder+1);
+    prefix << dir << OrderSuffix << (BondOrder+1);
     cout << "Writing files " << argv[1] << prefix.str() << ". ..." << endl;
     // energy
@@ -126 +142 @@
   // fragments
   prefix.str(" ");
-  prefix << argv[2] << EnergyFragmentSuffix;
+  prefix << dir << EnergyFragmentSuffix;
   if (!EnergyFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+  if (Hcorrected) {
+    prefix.str(" ");
+    prefix << dir << HcorrectionFragmentSuffix;
+    if (!HcorrectionFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
+  }
   prefix.str(" ");
-  prefix << argv[2] << ForceFragmentSuffix;
+  prefix << dir << ForceFragmentSuffix;
   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 (periode != NULL) { // also look for PAS values
     prefix.str(" ");
-    prefix << argv[2] << ShieldingFragmentSuffix;
+    prefix << dir << ShieldingFragmentSuffix;
     if (!ShieldingFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
     prefix.str(" ");
-    prefix << argv[2] << ShieldingPASFragmentSuffix;
+    prefix << dir << ShieldingPASFragmentSuffix;
     if (!ShieldingPASFragments.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], argv[2], EnergyFragmentSuffix)) return 1;
-  if (!Force.WriteLastMatrix(argv[1], argv[2], ForceFragmentSuffix)) return 1;
+  if (!Energy.WriteLastMatrix(argv[1], dir, EnergyFragmentSuffix)) return 1;
+  if (Hcorrected) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix);
+  if (!Force.WriteLastMatrix(argv[1], dir, ForceFragmentSuffix)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if (!Shielding.WriteLastMatrix(argv[1], argv[2], ShieldingFragmentSuffix)) return 1;
-    if (!ShieldingPAS.WriteLastMatrix(argv[1], argv[2], ShieldingPASFragmentSuffix)) return 1;
+    if (!Shielding.WriteLastMatrix(argv[1], dir, ShieldingFragmentSuffix)) return 1;
+    if (!ShieldingPAS.WriteLastMatrix(argv[1], dir, ShieldingPASFragmentSuffix)) return 1;
   }
 
   // exit  
   delete(periode);
+  Free((void **)&dir, "main: *dir");
   cout << "done." << endl;
   return 0;

src/moleculelist.cpp

@@ -133 +133 @@
 };
 
+/** Calculates necessary hydrogen correction due to unwanted interaction between saturated ones.
+ * If for a pair of two hydrogen atoms a and b, at least is a saturated one, and a and b are not
+ * bonded to the same atom, then we add for this pair a correction term constructed from a Morse
+ * potential function fit to QM calculations with respecting to the interatomic hydrogen distance.
+ * \param *out output stream for debugging
+ * \param *path path to file
+ */
+bool MoleculeListClass::AddHydrogenCorrection(ofstream *out, char *path)
+{
+  atom *Walker = NULL;
+  atom *Runner = NULL;
+  double ***FitConstant = NULL, **correction = NULL;
+  int a,b;
+  ofstream output;
+  ifstream input;
+  string line;
+  stringstream zeile;
+  double distance;
+  char ParsedLine[1023];
+  double tmp;
+  char *FragmentNumber = NULL;
+
+  cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
+  // 0. parse in fit constant files that should have the same dimension as the final energy files
+  // 0a. find dimension of matrices with constants
+  line = path;
+  line.append("/");
+  line += FRAGMENTPREFIX;
+  line += "1";
+  line += FITCONSTANTSUFFIX;
+  input.open(line.c_str());
+  if (input == NULL) {
+    cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
+    return false;
+  }
+  a=0;
+  b=-1; // we overcount by one 
+  while (!input.eof()) {
+    input.getline(ParsedLine, 1023);
+    zeile.str(ParsedLine);
+    int i=0;
+    while (!zeile.eof()) {
+      zeile >> distance;
+      i++;
+    }
+    if (i > a) 
+      a = i; 
+    b++;
+  }
+  cout << "I recognized " << a << " columns and " << b << " rows, ";
+  input.close();
+  
+  // 0b. allocate memory for constants
+  FitConstant = (double ***) Malloc(sizeof(double **)*3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  for (int k=0;k<3;k++) {
+    FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    for (int i=a;i--;) {
+      FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+    }
+  }
+  // 0c. parse in constants
+  for (int i=0;i<3;i++) {
+    line = path;
+    line.append("/");
+    line += FRAGMENTPREFIX;
+    sprintf(ParsedLine, "%d", i+1);
+    line += ParsedLine;
+    line += FITCONSTANTSUFFIX;
+    input.open(line.c_str());
+    if (input == NULL) {
+      cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
+      return false;
+    }
+    int k = 0,l;
+    while ((!input.eof()) && (k < b)) {
+      input.getline(ParsedLine, 1023);
+      //cout << "Current Line: " << ParsedLine << endl;
+      zeile.str(ParsedLine);
+      zeile.clear();
+      l = 0;
+      while ((!zeile.eof()) && (l < a)) {
+        zeile >> FitConstant[i][l][k];
+        //cout << FitConstant[i][l][k] << "\t";
+        l++;
+      }
+      //cout << endl;
+      k++;
+    }
+    input.close();
+  }
+  for(int k=0;k<3;k++) {
+    cout << "Constants " << k << ":" << endl;
+    for (int j=0;j<b;j++) {
+      for (int i=0;i<a;i++) {
+        cout << FitConstant[k][i][j] << "\t";
+      }
+      cout << endl;
+    }
+    cout << endl;
+  }
+  
+  // 0d. allocate final correction matrix
+  correction = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **correction");
+  for (int i=a;i--;)
+    correction[i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
+        
+  // 1a. go through every molecule in the list
+  for(int i=NumberOfMolecules;i--;) {
+    // 1b. zero final correction matrix
+    for (int k=a;k--;)
+      for (int j=b;j--;)
+        correction[k][j] = 0.;
+    // 2. take every hydrogen that is a saturated one
+    Walker = ListOfMolecules[i]->start;
+    while (Walker->next != ListOfMolecules[i]->end) {
+      Walker = Walker->next;
+      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
+      if ((Walker->type->Z == 1) && ((Walker->father == NULL) || (Walker->father->type->Z != 1))) { // if it's a hydrogen
+        Runner = ListOfMolecules[i]->start;
+        while (Runner->next != ListOfMolecules[i]->end) {
+          Runner = Runner->next;
+          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
+          // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
+          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) {  // (hydrogens have only one bonding partner!)
+            // 4. evaluate the morse potential for each matrix component and add up
+            distance = sqrt(Runner->x.Distance(&Walker->x));
+            //cout << "Fragment " << i << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
+            for(int k=0;k<a;k++) {
+              for (int j=0;j<b;j++) {
+                switch(k) {
+                case 1:
+                case 7:
+                case 11:
+                  tmp = pow(FitConstant[0][k][j] * ( 1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]) ) ),2);
+                  break;
+                default:
+                  tmp = FitConstant[0][k][j] * pow( distance,  FitConstant[1][k][j]) + FitConstant[2][k][j];
+                };
+                correction[k][j] -= tmp;    // ground state is actually lower (disturbed by additional interaction)
+                //cout << tmp << "\t";
+              }
+              //cout << endl;
+            }
+            //cout << endl;
+          }
+        }
+      }
+    }
+    // 5. write final matrix to file
+    line = path;
+    line.append("/");
+    line += FRAGMENTPREFIX;
+    FragmentNumber = FixedDigitNumber(NumberOfMolecules, i);
+    line += FragmentNumber;
+    delete(FragmentNumber);
+    line += HCORRECTIONSUFFIX;
+    output.open(line.c_str());
+    output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
+    for (int j=0;j<b;j++) {
+      for(int i=0;i<a;i++)
+        output << correction[i][j] << "\t";
+      output << endl;
+    }
+    output.close();
+  }
+  line = path;
+  line.append("/");
+  line += HCORRECTIONSUFFIX;
+  output.open(line.c_str());
+  output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
+  for (int j=0;j<b;j++) {
+    for(int i=0;i<a;i++)
+      output << 0 << "\t";
+    output << endl;
+  }
+  output.close();
+  // 6. free memory of parsed matrices
+  FitConstant = (double ***) Malloc(sizeof(double **)*a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  for (int k=0;k<3;k++) {
+    FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    for (int i=a;i--;) {
+      FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+    }
+  }
+  cout << "done." << endl;
+  return true;
+};
 
 /** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
@@ -224 +411 @@
     outputFragment.open(FragmentName, ios::out);
     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
-    if (intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment))
+    if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
       *out << " done." << endl;
     else
@@ -263 +450 @@
     configuration->SetDefaultPath(FragmentName);
     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
-    if (intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i]))
+    if ((intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
       *out << " done." << endl;
     else

src/molecules.cpp

@@ -2707 +2707 @@
     StoreAdjacencyToFile(out, configuration->configpath);
 
+    // store Hydrogen saturation correction file
+    BondFragments->AddHydrogenCorrection(out, configuration->configpath);
+    
     // store adaptive orders into file
     StoreOrderAtSiteFile(out, configuration->configpath);

src/molecules.hpp

@@ -289 +289 @@
 
   /// Output configs.
+  bool AddHydrogenCorrection(ofstream *out, char *path);
   bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
   bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);

src/parser.cpp

@@ -146 +146 @@
     file.str(" ");
     if (i != MatrixCounter)  
-      file << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix << suffix;
+      file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
     else
       file << name << prefix << suffix;
@@ -265 +265 @@
 };
 
+/** Scans all elements of MatrixContainer::Matrix for greatest absolute value.
+ * \return greatest value of MatrixContainer::Matrix
+ */
+double MatrixContainer::FindMaxValue()
+{
+  double max = Matrix[0][0][0];
+  for(int i=MatrixCounter+1;i--;)
+    for(int j=RowCounter[i]+1;j--;)
+      for(int k=ColumnCounter;k--;)
+        if (fabs(Matrix[i][j][k]) > max)
+          max = fabs(Matrix[i][j][k]);
+  if (fabs(max) < MYEPSILON)
+    max += MYEPSILON;
+ return max;
+};
+
+/** Scans all elements of MatrixContainer::Matrix for smallest absolute value.
+ * \return smallest value of MatrixContainer::Matrix
+ */
+double MatrixContainer::FindMinValue()
+{
+  double min = Matrix[0][0][0];
+  for(int i=MatrixCounter+1;i--;)
+    for(int j=RowCounter[i]+1;j--;)
+      for(int k=ColumnCounter;k--;)
+        if (fabs(Matrix[i][j][k]) < min)
+          min = fabs(Matrix[i][j][k]);
+  if (fabs(min) < MYEPSILON)
+    min += MYEPSILON;
+  return min;
+};
+
 /** Sets all values in the last of MatrixContainer::Matrix to \a value.
  * \param value reset value
@@ -333 +365 @@
               }
             }
-          //if ((ColumnCounter>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;
+            //if ((ColumnCounter>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 {
@@ -341 +373 @@
       }
     }
-    //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;
+   //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;
   }
   
@@ -428 +460 @@
  * Sums over the "F"-terms in ANOVA decomposition.
  * \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
  * \param Order bond order
+ * \parsm sign +1 or -1
  * \return true if summing was successful
  */
-bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, int Order)
+bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign)
 {
   // sum energy
-  for(int i=KeySet.FragmentsPerOrder[Order];i--;)
-    for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
-      for(int k=ColumnCounter;k--;)
-        Matrix[MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
+  if (CorrectionFragments == NULL)
+    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
+      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
+        for(int k=ColumnCounter;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--;)
+          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
   return true;
 };
@@ -492 +532 @@
  * \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 ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order)
-{
+bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
+{
+  int FragmentNr;
   // sum forces
-  for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++)
-    for(int l=0;l<RowCounter[ KeySet.OrderSet[Order][i] ];l++) {
-      int j = Indices[ KeySet.OrderSet[Order][i] ][l];
+  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 force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
         return false;
       }
-      if (j != -1)
+      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++)
-          Matrix[MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[Order][i] ][l][k];
-    }
+          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
+      }
+    }
+  }
   return true;
 };
@@ -571 +617 @@
       KeySets[i][j] = -1;
     FragmentNumber = FixedDigitNumber(FragmentCounter, i);
-    cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
+    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
     Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
     input.getline(filename, 1023);
@@ -577 +623 @@
     for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
       line >> KeySets[i][j];
-      cout << " " << KeySets[i][j];
-    }
-    cout << endl;
+      //cout << " " << KeySets[i][j];
+    }
+    //cout << endl;
   }
   input.close();

src/parser.hpp

@@ -19 +19 @@
 
 #define EnergySuffix ".energy.all"
+#define HcorrectionSuffix ".Hcorrection.all"
 #define ForcesSuffix ".forces.all"
 #define ShieldingSuffix ".sigma_all.csv"
@@ -26 +27 @@
 #define TimeSuffix ".speed"
 #define EnergyFragmentSuffix ".energyfragment.all"
+#define HcorrectionFragmentSuffix ".Hcorrectionfragment.all"
 #define ForceFragmentSuffix ".forcefragment.all"
 #define OrderSuffix ".Order"
@@ -52 +54 @@
   bool AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int CCounter);
   bool ResetMatrix();
+  double FindMinValue();
+  double FindMaxValue();
   bool SetLastMatrix(double value, int skipcolumns);
   bool SetLastMatrix(double **values, int skipcolumns);
@@ -66 +70 @@
   public:
     bool ParseIndices();
-    bool SumSubEnergy(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, int Order);
+    bool SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign);
 };
 
@@ -74 +78 @@
   public: 
     bool ParseIndices(char *name);
-    bool SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order);
+    bool SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign);
 };