Changes in src/joiner.cpp [437922:5bc4d0]

File:

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

src/joiner.cpp

@@ -2 +2 @@
  *
  * Takes evaluated fragments (energy and forces) and by reading the factors files determines total energy
- * and each force for the whole molecule.
- *
+ * and each force for the whole molecule. 
+ *    
  */
 
 //============================ INCLUDES ===========================
 
-#include "datacreator.hpp"
-#include "helpers.hpp"
-#include "parser.hpp"
-#include "periodentafel.hpp"
+#include "datacreator.hpp" 
+#include "helpers.hpp" 
+#include "parser.hpp" 
+#include "periodentafel.hpp" 
 
 //============================== MAIN =============================
@@ -28 +28 @@
   ForceMatrix ShieldingFragments;
   ForceMatrix ShieldingPASFragments;
-  KeySetsContainer KeySet;
+  EnergyMatrix Chi;
+  EnergyMatrix ChiPAS;
+  EnergyMatrix ChiFragments;
+  EnergyMatrix ChiPASFragments;
+  KeySetsContainer KeySet;  
   stringstream prefix;
   char *dir = NULL;
@@ -35 +39 @@
   cout << "Joiner" << endl;
   cout << "======" << endl;
-
+  
   // Get the command line options
   if (argc < 3) {
@@ -52 +56 @@
     periode->LoadPeriodentafel(argv[3]);
   }
-
+  
   // Test the given directory
   if (!TestParams(argc, argv))
     return 1;
-
+  
   // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++
-
+  
   // ------------- Parse through all Fragment subdirs --------
   if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix, 0,0)) return 1;
@@ -66 +70 @@
     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;
   }
 
@@ -71 +77 @@
   if (!Energy.ParseIndices()) return 1;
   if (Hcorrected) Hcorrection.ParseIndices();
-
+  
   // ---------- Parse the Force indices into an array ---------------
   if (!Force.ParseIndices(argv[1])) return 1;
@@ -79 +85 @@
     if(!Shielding.ParseIndices(argv[1])) return 1;
     if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
+    if(!Chi.ParseIndices()) return 1;
+    if(!ChiPAS.ParseIndices()) return 1;
   }
 
   // ---------- Parse the KeySets into an array ---------------
   if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1;
-
   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);
@@ -91 +99 @@
     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;
+  }
+  
   // ----------- Resetting last matrices (where full QM values are stored right now)
   if(!Energy.SetLastMatrix(0., 0)) return 1;
@@ -99 +109 @@
     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 +120 @@
     cout << "Summing energy of order " << BondOrder+1 << " ..." << endl;
     if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return 1;
-    if (Hcorrected) {
+    if (Hcorrected) { 
       HcorrectionFragments.SumSubManyBodyTerms(Hcorrection, KeySet, BondOrder);
       if (!Energy.SumSubEnergy(EnergyFragments, &HcorrectionFragments, KeySet, BondOrder, 1.)) return 1;
       if (Hcorrected) Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.);
-    } else
+    } else 
       if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1;
     // --------- sum up Forces --------------------
@@ -119 +131 @@
     if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return 1;
     if (periode != NULL) { // also look for PAS values
-      cout << "Summing shieldings of order " << BondOrder+1 << " ..." << endl;
+      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.SumSubEnergy(ChiFragments, NULL, KeySet, BondOrder, 1.)) return 1;
+      if (!ChiPASFragments.SumSubManyBodyTerms(ChiPAS, KeySet, BondOrder)) return 1;
+      if (!ChiPAS.SumSubEnergy(ChiPASFragments, NULL,KeySet, BondOrder, 1.)) return 1;
     }
 
@@ -138 +154 @@
       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;
     }
   }
@@ -160 +178 @@
     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;
   }
 
@@ -169 +193 @@
     if (!Shielding.WriteLastMatrix(argv[1], dir, ShieldingFragmentSuffix)) return 1;
     if (!ShieldingPAS.WriteLastMatrix(argv[1], dir, ShieldingPASFragmentSuffix)) return 1;
-  }
-
-  // exit
+    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");