Changes in molecuilder/src/analyzer.cpp [97c751:6ce722]

File:

• molecuilder/src/analyzer.cpp (modified) (14 diffs, 1 prop)

molecuilder/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;
-  ForceMatrix Chi;
-  ForceMatrix ChiPAS;
   EnergyMatrix Time;
-  EnergyMatrix EnergyFragments;
-  EnergyMatrix HcorrectionFragments;
-  ForceMatrix ForceFragments;
   ForceMatrix ShieldingFragments;
   ForceMatrix ShieldingPASFragments;
-  ForceMatrix ChiFragments;
-  ForceMatrix ChiPASFragments;
   KeySetsContainer KeySet;
   ofstream output;
@@ -53 +51 @@
   stringstream yrange;
   char *dir = NULL;
-  bool Hcorrected = true;
+  bool NoHCorrection = false;
+  bool NoHessian = false;
+  bool NoTime = false;
   double norm;
   int counter;
@@ -88 +88 @@
   // ------------- 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;
     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;
   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 ---------------
@@ -115 +132 @@
     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;
   }
 
@@ -129 +140 @@
   // ---------- 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;
     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;
   }
 
@@ -144 +156 @@
   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 +164 @@
   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 +172 @@
   output << endl;
 
-  if (periode != NULL) { // also look for PAS values
-    output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header << 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[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 +191 @@
     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;
     }
     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++) {
-    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 +223 @@
   // +++++++++++++++++++++++++++++++++++++++ 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 +284 @@
     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;
     }
-    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();
 
   
@@ -297 +316 @@
   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 +365 @@
   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
@@ -445 +466 @@
     }
     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();