Changes in / [cc2ee5:a98603]

Location:

src

Files:

• analyzer.cpp (modified) (8 diffs)
• joiner.cpp (modified) (9 diffs)
• parser.hpp (modified) (1 diff)

src/analyzer.cpp

@@ -29 +29 @@
   ForceMatrix Shielding;
   ForceMatrix ShieldingPAS;
+  ForceMatrix Chi;
+  ForceMatrix ChiPAS;
   EnergyMatrix Time;
   EnergyMatrix EnergyFragments;
@@ -35 +37 @@
   ForceMatrix ShieldingFragments;
   ForceMatrix ShieldingPASFragments;
+  ForceMatrix ChiFragments;
+  ForceMatrix ChiPASFragments;
   KeySetsContainer KeySet;
   ofstream output;
@@ -90 +94 @@
     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;
   }
 
@@ -109 +115 @@
     if(!ShieldingFragments.ParseIndices(argv[1])) return 1;
     if(!ShieldingPASFragments.ParseIndices(argv[1])) return 1;
+    if(!Chi.ParseIndices(argv[1])) return 1;
+    if(!ChiPAS.ParseIndices(argv[1])) return 1;
+    if (!ChiFragments.AllocateMatrix(Chi.Header, Chi.MatrixCounter, Chi.RowCounter, Chi.ColumnCounter)) return 1;
+    if (!ChiPASFragments.AllocateMatrix(ChiPAS.Header, ChiPAS.MatrixCounter, ChiPAS.RowCounter, ChiPAS.ColumnCounter)) return 1;
+    if(!ChiFragments.ParseIndices(argv[1])) return 1;
+    if(!ChiPASFragments.ParseIndices(argv[1])) return 1;
   }
 
@@ -122 +134 @@
     if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
     if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASFragmentSuffix, 1, 0)) return 1;
+    if (!ChiFragments.ParseFragmentMatrix(argv[1], dir, ChiFragmentSuffix, 1, 0)) return 1;
+    if (!ChiPASFragments.ParseFragmentMatrix(argv[1], dir, ChiPASFragmentSuffix, 1, 0)) return 1;
   }
 
@@ -159 +173 @@
       for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
         output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output << endl;
+
+    output << endl << "Total Chis" << endl << "===============" << endl << Chi.Header[Hessian.MatrixCounter] << endl;
+    for(int j=0;j<Chi.RowCounter[Chi.MatrixCounter];j++) {
+      for(int k=0;k<Chi.ColumnCounter[Chi.MatrixCounter];k++)
+        output << scientific << Chi.Matrix[ Chi.MatrixCounter ][j][k] << "\t";
+      output << endl;
+    }
+    output << endl;
+  
+    output << endl << "Total Chis PAS" << endl << "===============" << endl << ChiPAS.Header[ChiPAS.MatrixCounter] << endl;
+    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
+      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
+        output << scientific << ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t";
       output << endl;
     }
@@ -228 +258 @@
       output << endl;
     }
-  }
-  output.close();
+    output.close();
+    if (!CreateDataDeltaForcesOrderPerAtom(ChiPAS, ChiPASFragments, KeySet, argv[3], "DeltaChisPAS-Order", "Plot of error between approximated Chis and full Chis versus the Bond Order", datum)) return 1;
+    if (!CreateDataForcesOrderPerAtom(ChiPASFragments, KeySet, argv[3], "ChisPAS-Order", "Plot of approximated Chis versus the Bond Order", datum)) return 1;
+    if (!AppendOutputFile(output, argv[3], "ChisPAS-Order.dat" )) return false;
+    output << endl << "# Full" << endl;
+    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
+      output << j << "\t";
+      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
+        output << scientific <<  ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
+      output << endl;
+    }
+    output.close();
+  }
 
   
@@ -404 +445 @@
     }
     output << "'ShieldingsPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
+    output.close();  
+    output2.close();  
+
+    if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
+    if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
+    CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
+    CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
+    output << "set boxwidth " << step << endl;
+    output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
+    output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
+    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+      output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
+      output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
+      if (BondOrder-1 != KeySet.Order)
+        output2 << ", \\" << endl;
+    }
+    output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
     output.close();  
     output2.close();  

src/joiner.cpp

@@ -28 +28 @@
   ForceMatrix ShieldingFragments;
   ForceMatrix ShieldingPASFragments;
+  EnergyMatrix Chi;
+  EnergyMatrix ChiPAS;
+  EnergyMatrix ChiFragments;
+  EnergyMatrix ChiPASFragments;
   KeySetsContainer KeySet;  
   stringstream prefix;
@@ -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;
   }
 
@@ -79 +85 @@
     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;
   }
 
   // ---------- 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;
   }
   
@@ -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;
   }
 
@@ -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.SumSubForces(ChiFragments, KeySet, BondOrder, 1.)) return 1;
+      if (!ChiPASFragments.SumSubManyBodyTerms(ChiPAS, KeySet, BondOrder)) return 1;
+      if (!ChiPAS.SumSubForces(ChiPASFragments, 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;
+    if (!Chi.WriteLastMatrix(argv[1], dir, ChiFragmentSuffix)) return 1;
+    if (!ChiPAS.WriteLastMatrix(argv[1], dir, ChiPASFragmentSuffix)) return 1;
   }
 

src/parser.hpp

@@ -25 +25 @@
 #define ShieldingFragmentSuffix ".sigma_all_fragment.all"
 #define ShieldingPASFragmentSuffix ".sigma_all_PAS_fragment.all"
+#define ChiSuffix ".chi_all.csv"
+#define ChiPASSuffix ".chi_all_PAS.csv"
+#define ChiFragmentSuffix ".chi_all_fragment.all"
+#define ChiPASFragmentSuffix ".chi_all_PAS_fragment.all"
 #define TimeSuffix ".speed"
 #define EnergyFragmentSuffix ".energyfragment.all"