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

File:

• src/analyzer.cpp (modified) (25 diffs)

src/analyzer.cpp

@@ -3 +3 @@
  * Takes evaluated fragments (energy and forces) and does evaluation of how sensible the BOSSANOVA
  * approach was, e.g. in the decay of the many-body-contributions.
- *    
+ *
  */
 
@@ -9 +9 @@
 
 #include "datacreator.hpp"
-#include "helpers.hpp" 
+#include "helpers.hpp"
 #include "parser.hpp"
-#include "periodentafel.hpp" 
+#include "periodentafel.hpp"
 
 // include config.h
@@ -29 +29 @@
   ForceMatrix Shielding;
   ForceMatrix ShieldingPAS;
-  ForceMatrix Chi;
-  ForceMatrix ChiPAS;
   EnergyMatrix Time;
   EnergyMatrix EnergyFragments;
@@ -37 +35 @@
   ForceMatrix ShieldingFragments;
   ForceMatrix ShieldingPASFragments;
-  ForceMatrix ChiFragments;
-  ForceMatrix ChiPASFragments;
   KeySetsContainer KeySet;
   ofstream output;
@@ -54 +50 @@
   char *dir = NULL;
   bool Hcorrected = true;
-  double norm;
   int counter;
-  
+
   cout << "ANOVA Analyzer" << endl;
   cout << "==============" << endl;
-  
+
   // Get the command line options
   if (argc < 4) {
@@ -73 +68 @@
     strcat(dir, argv[2]);
   }
-  
+
   if (argc > 4) {
     cout << "Loading periodentafel." << endl;
@@ -79 +74 @@
     periode->LoadPeriodentafel(argv[4]);
   }
-  
+
   // Test the given directory
   if (!TestParams(argc, argv))
@@ -85 +80 @@
 
   // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++
-  
+
   // ------------- Parse through all Fragment subdirs --------
   if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
@@ -94 +89 @@
     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;
   }
 
@@ -101 +94 @@
   if (!Energy.ParseIndices()) return 1;
   if (Hcorrected) Hcorrection.ParseIndices();
-  
+
   // ---------- Parse the Force indices into an array ---------------
   if (!Force.ParseIndices(argv[1])) return 1;
@@ -115 +108 @@
     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;
   }
 
@@ -126 +113 @@
   if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1;
   if (!KeySet.ParseManyBodyTerms()) return 1;
-  
+
   // ---------- Parse fragment files created by 'joiner' into an array -------------
   if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
@@ -134 +121 @@
     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;
   }
 
   // +++++++++++++++ TESTING ++++++++++++++++++++++++++++++
-  
+
   // print energy and forces to file
   filename.str("");
@@ -151 +136 @@
   }
   output << endl;
-  
+
   output << endl << "Total Forces" << endl << "===============" << endl << Force.Header << endl;
   for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
@@ -168 +153 @@
     }
     output << endl;
-  
+
     output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header << endl;
     for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
@@ -176 +161 @@
     }
     output << endl;
-
-    output << endl << "Total Chis" << endl << "===============" << endl << Chi.Header[Chi.MatrixCounter] << endl;
-    for(int j=0;j<Chi.RowCounter[Chi.MatrixCounter];j++) {
-      for(int k=0;k<Chi.ColumnCounter;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;k++)
-        output << scientific << ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t";
-      output << endl;
-    }
-    output << endl;
-  }
-  
+  }
+
   output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
   for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
@@ -207 +176 @@
 
   // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
-  
+
   cout << "Analyzing ..." << endl;
 
@@ -258 +227 @@
       output << endl;
     }
-    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;k++)
-        output << scientific <<  ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
-      output << endl;
-    }
-    output.close();
-  }
-
-  
+  }
+  output.close();
+
+
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
   if (!CreateDataDeltaEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "DeltaEnergies-Order", "Plot of error between approximated and full energies energies versus the Bond Order", datum)) return 1;
-  
+
   // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
   if (!CreateDataEnergyOrder(EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
-  
+
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
   if (!CreateDataDeltaForcesOrderPerAtom(Force, ForceFragments, KeySet, argv[3], "DeltaForces-Order", "Plot of error between approximated forces and full forces versus the Bond Order", datum)) return 1;
@@ -341 +299 @@
 
   // ======================================= Creating the plot files ==============================================================
- 
+
   Orderxrange << "[1:" << KeySet.Order << "]";
   Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
   yrange.str("[1e-8:1e+1]");
-  
+
   // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
   if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "",  1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
-  
+
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
   if (!CreatePlotOrder(Energy, KeySet, argv[3], "DeltaEnergies-Order", 1, "outside", "y", "",  1, 1, "bond order k", "absolute error in energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
-  
+
   // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
   if (!CreatePlotOrder(Energy, KeySet, argv[3], "Energies-Order", 1, "outside", "", "",  1, 1, "bond order k", "approximate energy [Ht]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
@@ -369 +327 @@
   if(!OpenOutputFile(output, argv[3], "DeltaMinMeanMaxTotalForce-Order.pyx")) return 1;
   CreatePlotHeader(output, "DeltaMinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL,  1, 1, "bond order k", "absolute error in total forces [Ht/a.u.]");
-  output << "plot " << Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;  
+  output << "plot " << Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
   output << "'DeltaMinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
   output << "'DeltaMeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
   output << "'DeltaMaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
-  output.close();  
+  output.close();
 
   // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
@@ -388 +346 @@
   if(!OpenOutputFile(output, argv[3],"MinMeanMaxTotalForce-Order.pyx")) return 1;
   CreatePlotHeader(output, "MinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL, 1, 1, "bond order k", "absolute total force [Ht/a.u.]");
-  output << "plot "<< Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;  
+  output << "plot "<< Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;
   output << "'MinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
   output << "'MeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
   output << "'MaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
-  output.close();  
+  output.close();
 
   // ++++++++++++++++++++++++++++++++++++++Plotting vector sum vs. Order
@@ -411 +369 @@
   // min
   if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "minimum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
-    
+
   // mean
   if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "mean of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
-    
+
   // max
   if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "maximum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
@@ -421 +379 @@
   // min
   if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "minimum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
-    
+
   // mean
   if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "mean of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
-    
+
   // max
   if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "maximum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
- 
+
   // +++++++++++++++++++++++++++++++=Ploting approximated and true shielding for each atom
   if (periode != NULL) { // also look for PAS values
@@ -445 +403 @@
     }
     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();  
+    output.close();
+    output2.close();
   }
 
@@ -477 +418 @@
   output << "\trm -rf $(EPS)" << endl;
   output.close();
- 
+
   // ++++++++++++++++ exit ++++++++++++++++++++++++++++++++++
   delete(periode);