Changeset b38b64

Timestamp:

Jul 23, 2009, 9:14:18 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:

631dcb

Parents:

b77416 (diff), 72744a (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.

Message:

Merge branch 'HessianMatrix'

Location:

src

Files:

• analyzer.cpp (modified) (15 diffs)
• datacreator.cpp (modified) (36 diffs)
• datacreator.hpp (modified) (1 diff)
• joiner.cpp (modified) (12 diffs)
• parser.cpp (modified) (24 diffs)
• parser.hpp (modified) (5 diffs)

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;
@@ -32 +37 @@
   ForceMatrix ChiPAS;
   EnergyMatrix Time;
-  EnergyMatrix EnergyFragments;
-  EnergyMatrix HcorrectionFragments;
-  ForceMatrix ForceFragments;
   ForceMatrix ShieldingFragments;
   ForceMatrix ShieldingPASFragments;
@@ -53 +55 @@
   stringstream yrange;
   char *dir = NULL;
-  bool Hcorrected = true;
+  bool NoHCorrection = false;
+  bool NoHessian = false;
+  bool NoTime = false;
   double norm;
   int counter;
@@ -88 +92 @@
   // ------------- 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;
@@ -99 +114 @@
 
   // ---------- 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;
   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 ---------------
@@ -129 +152 @@
   // ---------- 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;
@@ -144 +170 @@
   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;
@@ -152 +178 @@
   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;
@@ -160 +186 @@
   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;
@@ -169 +205 @@
     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;
@@ -194 +230 @@
   }
   
-  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;
+  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 ++++++++++++++++++++++++++++++
@@ -214 +253 @@
   // +++++++++++++++++++++++++++++++++++++++ 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
@@ -254 +314 @@
     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;
@@ -297 +357 @@
   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;
@@ -346 +406 @@
   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

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.;
@@ -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, 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[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, char *dir, char *prefix, char *msg, 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, 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[ 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 +394 @@
   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 +418 @@
     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 +435 @@
     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 +441 @@
     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 +459 @@
   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 +465 @@
     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 +479 @@
 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 +490 @@
 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 +525 @@
 {
   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 +549 @@
 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 +581 @@
 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 +659 @@
 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 +683 @@
 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 +707 @@
 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 +715 @@
   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 +721 @@
     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 +738 @@
 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 +746 @@
   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 +752 @@
     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 +769 @@
 void BoxesForcePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
 {
-  stringstream line(Force.Header);
+  stringstream line(Force.Header[Force.MatrixCounter]);
   char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
   string token;
@@ -657 +778 @@
   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 +784 @@
     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 +801 @@
 void BoxesFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
 {
-  stringstream line(Force.Header);
+  stringstream line(Force.Header[Force.MatrixCounter]);
   char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
   string token;
@@ -689 +810 @@
   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 +816 @@
     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

@@ -26 +26 @@
 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 CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 
+bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 
+bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &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));

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;
@@ -35 +41 @@
   stringstream prefix;
   char *dir = NULL;
-  bool Hcorrected = true;
+  bool NoHCorrection = false;
+  bool NoHessian = false;
 
   cout << "Joiner" << endl;
@@ -65 +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;
@@ -75 +89 @@
 
   // ---------- 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 ---------------
@@ -94 +113 @@
 
   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;
@@ -106 +128 @@
   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;
@@ -120 +144 @@
     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.);
+      Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.);
     } else 
       if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1;
@@ -130 +154 @@
     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 and susceptibilities of order " << BondOrder+1 << " ..." << endl;
@@ -150 +180 @@
     // 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
@@ -162 +195 @@
   prefix << dir << EnergyFragmentSuffix;
   if (!EnergyFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
-  if (Hcorrected) {
+  if (!NoHCorrection) {
     prefix.str(" ");
     prefix << dir << HcorrectionFragmentSuffix;
@@ -171 +204 @@
   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(" ");
@@ -188 +226 @@
   // 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;

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][]");
@@ -89 +91 @@
   Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices"); 
   
-  Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
+  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.
@@ -120 +160 @@
     return false;
   }
-  
-  // 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
@@ -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
@@ -162 +203 @@
     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
+ * \param *CCounter number of columns for each matrix
  * \return Allocation successful 
  */
-bool MatrixContainer::AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int CCounter)
-{
-  Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader");
-  strncpy(Header, GivenHeader, 1023);
-
+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;
@@ -467 +513 @@
 // ======================================= CLASS EnergyMatrix =============================
 
-/** Create a trivial energy index mapping.
- * This just maps 1 to 1, 2 to 2 and so on for all fragments.
- * \return creation sucessful
- */
-bool EnergyMatrix::ParseIndices() 
-{
-  cout << "Parsing energy indices." << endl;
-  Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "EnergyMatrix::ParseIndices: **Indices");
-  for(int i=MatrixCounter+1;i--;) {
-    Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "EnergyMatrix::ParseIndices: *Indices[]");
-    for(int j=RowCounter[i];j--;)
-      Indices[i][j] = j;
-  } 
-  return true; 
-};
-
 /** 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 +528 @@
     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 +554 @@
     // 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
@@ -589 +623 @@
 
 /** 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 +643 @@
       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];
       }
@@ -655 +689 @@
     RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
     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(char *name, 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"
@@ -30 +36 @@
 #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 ==========================================
@@ -47 +49 @@
     double ***Matrix;
     int **Indices;
-    char *Header;
+    char **Header;
     int MatrixCounter;
     int *RowCounter;
-    int ColumnCounter;
+    int *ColumnCounter;
   
   MatrixContainer();
   ~MatrixContainer();
   
+  bool InitialiseIndices(class MatrixContainer *Matrix = NULL);
   bool ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr);
   virtual bool ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns);
-  bool AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int CCounter);
+  bool AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter);
   bool ResetMatrix();
   double FindMinValue();
@@ -74 +77 @@
 class EnergyMatrix : public MatrixContainer {
   public:
-    bool ParseIndices();
     bool SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign);
     bool ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns);
@@ -86 +88 @@
     bool SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign);
     bool ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns);
+};
+
+// ======================================= 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(char *name, char *prefix, string suffix, int skiplines, int skipcolumns);
+  private:
+    bool IsSymmetric;
 };