Ignore:
Timestamp:
Jul 24, 2008, 2:00:19 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
747b10
Parents:
b3eb8e
Message:

BIG change: Hcorrection included and bugfix of analyzer (wrt forces)

  • Hcorrection is supplied via molecuilder which gives correction values if hydrogen atoms are to close and thus already feel each other's influence. This correction energy matrix is also parsed ans subtracted from the approximated energy matrix generated from the fragments. This impacts both on joiner, analyzer and molecuilder
  • Forces were not working which had multiple causes:
    • first, we stepped back down to absolute errors which rather gives a feel about convergence (1e-6 is simply small and neglegible)
    • there may have been bugs where (either in Force or in ForceFragments) adding up took place, this is cleaned up
    • more routines have been abstracted from analyzer into datacreator functions
      • the order of ions had changed, re-run with current molecuilder, calculation and subsequent joining/analyzing brought correct results finally
  • plots were bad still due to wrong axes (two more functions finding max/min values), use of wrong columns in plot file
File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/analyzer.cpp

    rb3eb8e rdb3ea3  
    2525  periodentafel *periode = NULL; // and a period table of all elements
    2626  EnergyMatrix Energy;
     27  EnergyMatrix Hcorrection;
    2728  ForceMatrix Force;
    2829  ForceMatrix Shielding;
    … …  
    3031  MatrixContainer Time;
    3132  EnergyMatrix EnergyFragments;
     33  EnergyMatrix HcorrectionFragments;
    3234  ForceMatrix ForceFragments;
    3335  ForceMatrix ShieldingFragments;
    … …  
    4547  stringstream Orderxrange;
    4648  stringstream Fragmentxrange;
    47 
     49  stringstream yrange;
     50  char *dir = NULL;
     51  bool Hcorrected = true;
     52  double norm;
     53 
    4854  cout << "ANOVA Analyzer" << endl;
    4955  cout << "==============" << endl;
    … …  
    5763    cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl;
    5864    return 1;
    59   }
     65  } else {
     66    dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir");
     67    strcpy(dir, "/");
     68    strcat(dir, argv[2]);
     69  }
     70 
    6071  if (argc > 4) {
    6172    cout << "Loading periodentafel." << endl;
    … …  
    7182 
    7283  // ------------- Parse through all Fragment subdirs --------
    73   if (!Energy.ParseMatrix(argv[1], argv[2], EnergySuffix,0,0)) return 1;
    74   if (!Force.ParseMatrix(argv[1], argv[2], ForcesSuffix,0,0)) return 1;
    75   if (!Time.ParseMatrix(argv[1], argv[2], TimeSuffix, 10,1)) return 1;
     84  if (!Energy.ParseMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
     85  Hcorrected = Hcorrection.ParseMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
     86  if (!Force.ParseMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
     87  //if (!Time.ParseMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
    7688  if (periode != NULL) { // also look for PAS values
    77     if (!Shielding.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
    78     if (!ShieldingPAS.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
     89    if (!Shielding.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
     90    if (!ShieldingPAS.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
    7991  }
    8092
    8193  // ---------- Parse the TE Factors into an array -----------------
    8294  if (!Energy.ParseIndices()) return 1;
     95  if (Hcorrected) Hcorrection.ParseIndices();
    8396 
    8497  // ---------- Parse the Force indices into an array ---------------
    8598  if (!Force.ParseIndices(argv[1])) return 1;
     99  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
     100  if (!ForceFragments.ParseIndices(argv[1])) return 1;
    86101
    87102  // ---------- Parse the shielding indices into an array ---------------
    … …  
    96111 
    97112  // ---------- Parse fragment files created by 'joiner' into an array -------------
    98   if (!EnergyFragments.ParseMatrix(argv[1], argv[2], EnergyFragmentSuffix,0,0)) return 1;
    99   if (!ForceFragments.ParseMatrix(argv[1], argv[2], ForceFragmentSuffix,0,0)) return 1;
     113  if (!EnergyFragments.ParseMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
     114  if (Hcorrected) HcorrectionFragments.ParseMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
     115  if (!ForceFragments.ParseMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
    100116  if (periode != NULL) { // also look for PAS values
    101     if (!ShieldingFragments.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
    102     if (!ShieldingPASFragments.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
     117    if (!ShieldingFragments.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
     118    if (!ShieldingPASFragments.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
    103119  }
    104120
    … …  
    143159  }
    144160 
    145   output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
    146   Time.SetLastMatrix(0., 0);
    147   for (int BondOrder=KeySet.Order;BondOrder--;)
    148     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
    149       for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
    150         for(int k=Time.ColumnCounter;k--;) {
    151           Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    152         }
    153   for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
    154     for(int k=0;k<Time.ColumnCounter;k++)
    155       output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
    156     output << endl;
    157   }
    158   output << endl;
     161//  output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
     162//  Time.SetLastMatrix(0., 0);
     163//  for (int BondOrder=KeySet.Order;BondOrder--;)
     164//    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
     165//      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
     166//        for(int k=Time.ColumnCounter;k--;) {
     167//          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     168//        }
     169//  for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
     170//    for(int k=0;k<Time.ColumnCounter;k++)
     171//      output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
     172//    output << endl;
     173//  }
     174//  output << endl;
    159175  output.close();
    160176
    … …  
    165181  // ======================================= Creating the data files ==============================================================
    166182
    167   // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    168   if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
    169   Time.SetLastMatrix(0., 0);
    170   output << "#Order\tFrag.No.\t" << Time.Header << endl;
    171   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    172     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
    173       for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
    174         for(int k=Time.ColumnCounter;k--;) {
    175           Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    176         }
    177     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    178     for(int k=0;k<Time.ColumnCounter;k++)
    179       output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][1][k];
    180     output << endl;
    181   }
    182   output.close();
     183//  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
     184//  if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
     185//  Time.SetLastMatrix(0., 0);
     186//  output << "#Order\tFrag.No.\t" << Time.Header << endl;
     187//  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     188//    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
     189//      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
     190//        for(int k=Time.ColumnCounter;k--;) {
     191//          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     192//        }
     193//    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     194//    for(int k=0;k<Time.ColumnCounter;k++)
     195//      output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][1][k];
     196//    output << endl;
     197//  }
     198//  output.close();
    183199
    184200  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
    185   if (!EnergyFragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0)) return 1;
    186   if (!CreateDataEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "DeltaEnergies-Order", "Plot of error between approximated and full energies energies versus the Bond Order", datum)) return 1;
     201  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;
    187202 
    188203  // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
    189   if (!EnergyFragments.SetLastMatrix(1.,0)) return 1;
    190   if (!CreateDataEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
     204  if (!CreateDataEnergyOrder(EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
    191205 
    192206  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
    193   if (!ForceFragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0)) return 1;
    194   if (!OpenOutputFile(output, argv[3], "DeltaForces-Order.dat" )) return false;
    195   if (!OpenOutputFile(output2, argv[3], "DeltaMinForces-Order.dat" )) return false;
    196   if (!OpenOutputFile(output3, argv[3], "DeltaMeanForces-Order.dat" )) return false;
    197   if (!OpenOutputFile(output4, argv[3], "DeltaMaxForces-Order.dat" )) return false;
    198   cout << "Plot of per atom and min/mean/max error between approximated forces and full forces versus the Bond Order ... " << endl;
    199   output << "# Plot of error between approximated forces and full forces versus the Bond Order, created on " << datum;
    200   output << "# AtomNo" << Force.Header << endl;
    201   output2 << "# Plot of min error between approximated forces and full forces versus the Bond Order, created on " << datum;
    202   output2 << "# Order\tFrag.No.\t" << Force.Header << endl;
    203   output3 << "# Plot of mean error between approximated forces and full forces versus the Bond Order, created on " << datum;
    204   output3 << "# Order\tFrag.No.\t" << Force.Header << endl;
    205   output4 << "# Plot of max error between approximated forces and full forces versus the Bond Order, created on " << datum;
    206   output4 << "# Order\tFrag.No.\t" << Force.Header << endl;
    207   int FragmentCounter = 0;
    208   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    209     for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
    210       for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
    211         int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
    212         if (j > Force.RowCounter[Force.MatrixCounter]) {
    213           cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
    214           return 1;
    215         }
    216         if (j != -1)
    217           for(int k=0;k<Force.ColumnCounter;k++) {
    218             Force.Matrix[Force.MatrixCounter][j][k] -= ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
    219           }
    220       }
    221       FragmentCounter++;
    222     }
    223     // errors per atom
    224     output << "#Order\t" << BondOrder+1 << endl;
    225     for(int j=0;j<Force.RowCounter[ Force.MatrixCounter ];j++) {
    226       output << Force.Indices[Force.MatrixCounter][j] << "\t";
    227       for (int l=0;l<Force.ColumnCounter;l++)
    228         if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ j ][l]) < MYEPSILON)
    229           output << scientific << Force.Matrix[Force.MatrixCounter][ j ][l] << "\t";
    230         else
    231           output << scientific << (Force.Matrix[Force.MatrixCounter][ j ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ j ][l]) << "\t";
    232       output << endl;
    233     }
    234     output << endl;
    235    
    236     // min error
    237     output2 << BondOrder+1 << "\t" << FragmentCounter;
    238     CreateMinimumForce(Force, Force.MatrixCounter);
    239     for (int l=0;l<Force.ColumnCounter;l++)
    240       if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
    241         output2 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
    242       else
    243         output2 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
    244     output2 << endl;
    245    
    246     // mean error
    247     output3 << BondOrder+1 << "\t" << FragmentCounter;
    248     CreateMeanForce(Force, Force.MatrixCounter);
    249     for (int l=0;l<Force.ColumnCounter;l++)
    250       if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
    251         output3 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
    252       else
    253         output3 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
    254     output3 << endl;
    255    
    256     // max error
    257     output4 << BondOrder+1 << "\t" << FragmentCounter;
    258     CreateMaximumForce(Force, Force.MatrixCounter);
    259     for (int l=0;l<Force.ColumnCounter;l++)
    260       if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
    261         output4 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
    262       else
    263         output4 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
    264     output4 << endl;
    265   }
    266   output.close();
    267   output2.close();
    268   output3.close();
    269   output4.close();
     207  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;
     208
     209  // min force
     210  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMinForces-Order", "Plot of min error between approximated forces and full forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
     211
     212  // mean force
     213  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMeanForces-Order", "Plot of mean error between approximated forces and full forces versus the Bond Order", datum, CreateMeanForce)) return 1;
     214
     215  // max force
     216  if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMaxForces-Order", "Plot of max error between approximated forces and full forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
    270217
    271218  // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
    272   if(!OpenOutputFile(output, argv[3], "Forces-Order.dat")) return 1;
    273   cout << "Plot of approximated forces versus the Bond Order ... " << endl;
    274   output << "# Plot of approximated forces versus the Bond Order, created on " << datum;
    275   output << "# AtomNo" << Force.Header << endl;
    276   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    277     for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
    278       for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
    279         int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
    280         if (j > Force.RowCounter[Force.MatrixCounter]) {
    281           cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
    282           return 1;
    283         }
    284         if (j != -1)
    285           for(int k=0;k<Force.ColumnCounter;k++) {
    286             Force.Matrix[Force.MatrixCounter][j][k] -= ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
    287           }
    288       }
    289     }
    290     // errors per atom
    291     output << "#Order\t" << BondOrder+1 << endl;
    292     for(int j=0;j<Force.RowCounter[ Force.MatrixCounter ];j++) {
    293       output << Force.Indices[Force.MatrixCounter][j] << "\t";
    294       for (int l=0;l<Force.ColumnCounter;l++)
    295          output << Force.Matrix[Force.MatrixCounter][ j ][l] << "\t";
    296       output << endl;
    297     }
    298     output << endl;
    299   }
    300   output.close();
    301 
    302   // min force
    303   if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
    304 
    305   // mean force
    306   if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;
    307 
    308   // max force
    309   if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
     219  if (!CreateDataForcesOrderPerAtom(ForceFragments, KeySet, argv[3], "Forces-Order", "Plot of approximated forces versus the Bond Order", datum)) return 1;
     220
     221  // min force
     222  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;
     223
     224  // mean force
     225  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;
     226
     227  // max force
     228  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;
    310229
    311230  // ++++++++++++++++++++++++++++++++++++++Plotting vector sum (should be 0) vs. bond order
    312   if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;
     231  if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;
    313232
    314233  // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
    … …  
    342261  Orderxrange << "[1:" << KeySet.Order << "]";
    343262  Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
     263  yrange.str("[1e-8:1e+1]");
    344264 
    345265  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    346   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;
     266  //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;
    347267 
    348268  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
    349   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(), "[1e-8:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
     269  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;
    350270 
    351271  // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
    … …  
    353273
    354274  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
    355   // min force
    356   if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMinForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in min force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
    357 
    358   // mean force
    359   if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMeanForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in mean force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;
    360 
    361   // max force
    362   if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMaxForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in max force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
     275  yrange.str("[1e-8:1e+0]");
     276  // min force
     277  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMinForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in min force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
     278
     279  // mean force
     280  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMeanForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in mean force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;
     281
     282  // max force
     283  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMaxForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in max force [Ht/a.u.]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;
    363284
    364285  // min/mean/max comparison for total force
    … …  
    395316
    396317  // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
    397   if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-Fragment", 5, "below", "y", "", 1, 5, "fragment number", "Energies of each fragment [Ht]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with points", AbsEnergyPlotLine)) return 1;
    398   if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Energies of each fragment [Ht]", Orderxrange.str().c_str(), "[1e-10:1e+1]", "1" , "with points", AbsEnergyPlotLine)) return 1;
    399   if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Maximum fragment energy [Ht]", Orderxrange.str().c_str(), "[1e-6:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
    400   if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Minimum fragment energy [Ht]", Orderxrange.str().c_str(), "[1e-6:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
     318  yrange.str("");
     319  yrange << "[" << EnergyFragments.FindMinValue() << ":" << EnergyFragments.FindMaxValue() << "]";
     320  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-Fragment", 5, "below", "y", "", 1, 5, "fragment number", "Energies of each fragment [Ht]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "2" , "with points", AbsEnergyPlotLine)) return 1;
     321  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Energies of each fragment [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with points", AbsEnergyPlotLine)) return 1;
     322  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Maximum fragment energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
     323  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Minimum fragment energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
    401324
    402325  // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment
     326  yrange.str("");
     327  yrange << "[" << ForceFragments.FindMinValue() << ":" << ForceFragments.FindMaxValue()<< "]";
    403328  // min
    404   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(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
     329  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;
    405330   
    406331  // mean
    407   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(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
     332  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;
    408333   
    409334  // max
    410   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(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
     335  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;
    411336
    412337  // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment per bond order
    413338  // min
    414   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(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
     339  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;
    415340   
    416341  // mean
    417   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(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
     342  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;
    418343   
    419344  // max
    420   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(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
     345  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;
    421346 
    422347  // create Makefile
    … …  
    434359  // ++++++++++++++++ exit ++++++++++++++++++++++++++++++++++
    435360  delete(periode);
     361  Free((void **)&dir, "main: *dir");
    436362  cout << "done." << endl;
    437363  return 0;
Note: See TracChangeset for help on using the changeset viewer.