Changes in src/joiner.cpp [042f82:29812d]

File:

• src/joiner.cpp (modified) (12 diffs)

src/joiner.cpp

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