Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/datacreator.cpp

    r6ac7ee rce5ac3  
    11/** \file datacreator.cpp
    22 *
    3  * Declarations of assisting functions in creating data and plot files.
    4  *
     3 * Declarations of assisting functions in creating data and plot files. 
     4 *   
    55 */
    66
    … …  
    1919bool OpenOutputFile(ofstream &output, const char *dir, const char *filename)
    2020{
    21         stringstream name;
    22         name << dir << "/" << filename;
    23         output.open(name.str().c_str(), ios::out);
    24         if (output == NULL) {
    25                 cout << "Unable to open " << name.str() << " for writing, is directory correct?" << endl;
    26                 return false;
    27         }
    28         return true;
    29 };
     21  stringstream name;
     22  name << dir << "/" << filename;
     23  output.open(name.str().c_str(), ios::out);
     24  if (output == NULL) {
     25    cout << "Unable to open " << name.str() << " for writing, is directory correct?" << endl;
     26    return false;
     27  }
     28  return true;
     29}; 
    3030
    3131/** Opens a file for appending with \a *filename in \a *dir.
    … …  
    3737bool AppendOutputFile(ofstream &output, const char *dir, const char *filename)
    3838{
    39         stringstream name;
    40         name << dir << "/" << filename;
    41         output.open(name.str().c_str(), ios::app);
    42         if (output == NULL) {
    43                 cout << "Unable to open " << name.str() << " for writing, is directory correct?" << endl;
    44                 return false;
    45         }
    46         return true;
    47 };
     39  stringstream name;
     40  name << dir << "/" << filename;
     41  output.open(name.str().c_str(), ios::app);
     42  if (output == NULL) {
     43    cout << "Unable to open " << name.str() << " for writing, is directory correct?" << endl;
     44    return false;
     45  }
     46  return true;
     47}; 
    4848
    4949/** Plots an energy vs. order.
    … …  
    5454 * \return true if file was written successfully
    5555 */
    56 bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
    57 {
    58         stringstream filename;
    59         ofstream output;
    60 
    61         filename << prefix << ".dat";
    62         if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
    63         cout << msg << endl;
    64         output << "# " << msg << ", created on " << datum;
    65         output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
    66         for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    67                 for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
    68                         for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
    69                                 for(int k=Fragments.ColumnCounter;k--;)
    70                                         Fragments.Matrix[Fragments.MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    71                 }
    72                 output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    73                 for (int l=0;l<Fragments.ColumnCounter;l++)
    74                         output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
    75                 output << endl;
    76         }
    77         output.close();
    78         return true;
     56bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
     57{
     58  stringstream filename;
     59  ofstream output;
     60
     61  filename << prefix << ".dat";
     62  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     63  cout << msg << endl;
     64  output << "# " << msg << ", created on " << datum;
     65  output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
     66  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     67    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
     68      for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
     69        for(int k=Fragments.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];k--;)
     70          Fragments.Matrix[Fragments.MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     71    }
     72    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     73    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
     74      output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
     75    output << endl;
     76  }
     77  output.close();
     78  return true;
    7979};
    8080
    … …  
    8787 * \return true if file was written successfully
    8888 */
    89 bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
    90 {
    91         stringstream filename;
    92         ofstream output;
    93 
    94         filename << prefix << ".dat";
    95         if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
    96         cout << msg << endl;
    97         output << "# " << msg << ", created on " << datum;
    98         output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
    99         Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0);
    100         for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    101                 for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
    102                         for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
    103                                 for(int k=Fragments.ColumnCounter;k--;)
    104                                         Fragments.Matrix[Fragments.MatrixCounter][j][k] -= Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    105                 }
    106                 output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    107                 for (int l=0;l<Fragments.ColumnCounter;l++)
    108                         if (fabs(Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]) < MYEPSILON)
    109                                 output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
    110                         else
    111                                 output << scientific << "\t" << (Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l] / Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]);
    112                 output << endl;
    113         }
    114         output.close();
    115         return true;
     89bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
     90{
     91  stringstream filename;
     92  ofstream output;
     93
     94  filename << prefix << ".dat";
     95  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     96  cout << msg << endl;
     97  output << "# " << msg << ", created on " << datum;
     98  output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
     99  Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0);
     100  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     101    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
     102      for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
     103        for(int k=Fragments.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];k--;)
     104          Fragments.Matrix[Fragments.MatrixCounter][j][k] -= Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     105    }
     106    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     107    for (int l=0;l<Fragments.ColumnCounter[Energy.MatrixCounter];l++)
     108      if (fabs(Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]) < MYEPSILON)
     109        output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
     110      else
     111        output << scientific << "\t" << (Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l] / Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]);
     112    output << endl;
     113  }
     114  output.close();
     115  return true;
    116116};
    117117
    … …  
    126126bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
    127127{
    128         stringstream filename;
    129         ofstream output;
    130 
    131         filename << prefix << ".dat";
    132         if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
    133         cout << msg << endl;
    134         output << "# " << msg << ", created on " << datum;
    135         output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
    136         Fragments.SetLastMatrix(0.,0);
    137         for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    138                 Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
    139                 output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    140                 CreateForce(Fragments, Fragments.MatrixCounter);
    141                 for (int l=0;l<Fragments.ColumnCounter;l++)
    142                         output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
    143                 output << endl;
    144         }
    145         output.close();
    146         return true;
     128  stringstream filename;
     129  ofstream output;
     130
     131  filename << prefix << ".dat";
     132  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     133  cout << msg << endl;
     134  output << "# " << msg << ", created on " << datum;
     135  output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
     136  Fragments.SetLastMatrix(0.,0);
     137  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     138    Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
     139    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     140    CreateForce(Fragments, Fragments.MatrixCounter);
     141    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
     142      output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
     143    output << endl;
     144  }
     145  output.close();
     146  return true;
    147147};
    148148
    … …  
    158158bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
    159159{
    160         stringstream filename;
    161         ofstream output;
    162 
    163         filename << prefix << ".dat";
    164         if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
    165         cout << msg << endl;
    166         output << "# " << msg << ", created on " << datum;
    167         output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
    168         Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
    169         for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    170                 Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
    171                 output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    172                 CreateForce(Fragments, Fragments.MatrixCounter);
    173                 for (int l=0;l<Fragments.ColumnCounter;l++)
    174                         output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
    175                 output << endl;
    176         }
    177         output.close();
    178         return true;
     160  stringstream filename;
     161  ofstream output;
     162
     163  filename << prefix << ".dat";
     164  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     165  cout << msg << endl;
     166  output << "# " << msg << ", created on " << datum;
     167  output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
     168  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
     169  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     170    Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
     171    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     172    CreateForce(Fragments, Fragments.MatrixCounter);
     173    for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
     174      output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
     175    output << endl;
     176  }
     177  output.close();
     178  return true;
    179179};
    180180
    … …  
    188188 * \return true if file was written successfully
    189189 */
    190 bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
    191 {
    192         stringstream filename;
    193         ofstream output;
    194         double norm = 0.;
    195 
    196         filename << prefix << ".dat";
    197         if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
    198         cout << msg << endl;
    199         output << "# " << msg << ", created on " << datum;
    200         output << "# AtomNo\t" << Fragments.Header << endl;
    201         Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
    202         for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    203                 //cout << "Current order is " << BondOrder << "." << endl;
    204                 Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
    205                 // errors per atom
    206                 output << endl << "#Order\t" << BondOrder+1 << endl;
    207                 for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
    208                         output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
    209                         for (int l=0;l<Fragments.ColumnCounter;l++) {
    210                                 if (((l+1) % 3) == 0) {
    211                                         norm = 0.;
    212                                         for (int m=0;m<NDIM;m++)
    213                                                 norm += Force.Matrix[Force.MatrixCounter][ j ][l+m]*Force.Matrix[Force.MatrixCounter][ j ][l+m];
    214                                         norm = sqrt(norm);
    215                                 }                                                                                                                                                                                                                                                                                                                                                                                                               
    216 //                              if (norm < MYEPSILON)
    217                                         output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
    218 //                              else
    219 //                                      output << scientific << (Fragments.Matrix[Fragments.MatrixCounter][ j ][l] / norm) << "\t";
    220                         }
    221                         output << endl;
    222                 }
    223                 output << endl;
    224         }
    225         output.close();
    226         return true;
     190bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 
     191{
     192  stringstream filename;
     193  ofstream output;
     194  double norm = 0.;
     195
     196  filename << prefix << ".dat";
     197  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     198  cout << msg << endl;
     199  output << "# " << msg << ", created on " << datum;
     200  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
     201  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
     202  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     203    //cout << "Current order is " << BondOrder << "." << endl;
     204    Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
     205    // errors per atom
     206    output << endl << "#Order\t" << BondOrder+1 << endl;
     207    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
     208      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
     209      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
     210        if (((l+1) % 3) == 0) {
     211          norm = 0.;
     212          for (int m=0;m<NDIM;m++)
     213            norm += Force.Matrix[Force.MatrixCounter][ j ][l+m]*Force.Matrix[Force.MatrixCounter][ j ][l+m];
     214          norm = sqrt(norm);
     215        }                                                                                                           
     216//        if (norm < MYEPSILON)
     217          output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
     218//        else
     219//          output << scientific << (Fragments.Matrix[Fragments.MatrixCounter][ j ][l] / norm) << "\t";
     220      }
     221      output << endl;
     222    }
     223    output << endl;
     224  }
     225  output.close();
     226  return true;
    227227};
    228228
    … …  
    235235 * \return true if file was written successfully
    236236 */
    237 bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
    238 {
    239         stringstream filename;
    240         ofstream output;
    241 
    242         filename << prefix << ".dat";
    243         if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
    244         cout << msg << endl;
    245         output << "# " << msg << ", created on " << datum;
    246         output << "# AtomNo\t" << Fragments.Header << endl;
    247         for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    248                 //cout << "Current order is " << BondOrder << "." << endl;
    249                 Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
    250                 // errors per atom
    251                 output << endl << "#Order\t" << BondOrder+1 << endl;
    252                 for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
    253                         output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
    254                         for (int l=0;l<Fragments.ColumnCounter;l++)
    255                                 output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
    256                         output << endl;
    257                 }
    258                 output << endl;
    259         }
    260         output.close();
    261         return true;
     237bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
     238{
     239  stringstream filename;
     240  ofstream output;
     241
     242  filename << prefix << ".dat";
     243  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     244  cout << msg << endl;
     245  output << "# " << msg << ", created on " << datum;
     246  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
     247  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     248    //cout << "Current order is " << BondOrder << "." << endl;
     249    Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
     250    // errors per atom
     251    output << endl << "#Order\t" << BondOrder+1 << endl;
     252    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
     253      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
     254      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++)
     255        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
     256      output << endl;
     257    }
     258    output << endl;
     259  }
     260  output.close();
     261  return true;
     262};
     263
     264
     265/** Plot hessian error vs. vs atom vs. order.
     266 * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
     267 * \param &Fragments HessianMatrix class containing matrix values
     268 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
     269 * \param *prefix prefix in filename (without ending)
     270 * \param *msg message to be place in first line as a comment
     271 * \param *datum current date and time
     272 * \return true if file was written successfully
     273 */
     274bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
     275{
     276  stringstream filename;
     277  ofstream output;
     278
     279  filename << prefix << ".dat";
     280  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     281  cout << msg << endl;
     282  output << "# " << msg << ", created on " << datum;
     283  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
     284  Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
     285  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     286    //cout << "Current order is " << BondOrder << "." << endl;
     287    Fragments.SumSubHessians(Fragments, KeySet, BondOrder, -1.);
     288    // errors per atom
     289    output << endl << "#Order\t" << BondOrder+1 << endl;
     290    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
     291      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
     292      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
     293        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
     294      }
     295      output << endl;
     296    }
     297    output << endl;
     298  }
     299  output.close();
     300  return true;
     301};
     302
     303/** Plot hessian error vs. vs atom vs. order in the frobenius norm.
     304 * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
     305 * \param &Fragments HessianMatrix class containing matrix values
     306 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
     307 * \param *prefix prefix in filename (without ending)
     308 * \param *msg message to be place in first line as a comment
     309 * \param *datum current date and time
     310 * \return true if file was written successfully
     311 */
     312bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
     313{
     314  stringstream filename;
     315  ofstream output;
     316  double norm = 0;
     317  double tmp;
     318
     319  filename << prefix << ".dat";
     320  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     321  cout << msg << endl;
     322  output << "# " << msg << ", created on " << datum;
     323  output << "# AtomNo\t";
     324  Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
     325  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     326    output << "Order" << BondOrder+1 << "\t";
     327  }
     328  output << endl;
     329  output << Fragments.RowCounter[ Fragments.MatrixCounter ] << "\t";
     330  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     331    //cout << "Current order is " << BondOrder << "." << endl;
     332    Fragments.SumSubHessians(Fragments, KeySet, BondOrder, -1.);
     333    // frobenius norm of errors per atom
     334    norm = 0.;
     335    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
     336      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
     337        tmp = Fragments.Matrix[Fragments.MatrixCounter][ j ][l];
     338        norm += tmp*tmp;
     339      }
     340    }
     341    output << scientific << sqrt(norm)/(Fragments.RowCounter[ Fragments.MatrixCounter ]*Fragments.ColumnCounter[ Fragments.MatrixCounter] ) << "\t";
     342  }
     343  output << endl;
     344  output.close();
     345  return true;
     346};
     347
     348/** Plot hessian error vs. vs atom vs. order.
     349 * \param &Fragments HessianMatrix class containing matrix values
     350 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
     351 * \param *prefix prefix in filename (without ending)
     352 * \param *msg message to be place in first line as a comment
     353 * \param *datum current date and time
     354 * \return true if file was written successfully
     355 */
     356bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
     357{
     358  stringstream filename;
     359  ofstream output;
     360
     361  filename << prefix << ".dat";
     362  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     363  cout << msg << endl;
     364  output << "# " << msg << ", created on " << datum;
     365  output << "# AtomNo\t" << Fragments.Header[ Fragments.MatrixCounter ] << endl;
     366  Fragments.SetLastMatrix(0., 0);
     367  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     368    //cout << "Current order is " << BondOrder << "." << endl;
     369    Fragments.SumSubHessians(Fragments, KeySet, BondOrder, 1.);
     370    // errors per atom
     371    output << endl << "#Order\t" << BondOrder+1 << endl;
     372    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
     373      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
     374      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++)
     375        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
     376      output << endl;
     377    }
     378    output << endl;
     379  }
     380  output.close();
     381  return true;
    262382};
    263383
    … …  
    266386bool CreateDataFragment(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateFragment)(class MatrixContainer &, int))
    267387{
    268         stringstream filename;
    269         ofstream output;
    270 
    271         filename << prefix << ".dat";
    272         if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
    273         cout << msg << endl;
    274         output << "# " << msg << ", created on " << datum << endl;
    275         output << "#Order\tFrag.No.\t" << Fragment.Header << endl;
    276         for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    277                 for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
    278                         output << BondOrder+1 << "\t" << KeySet.OrderSet[BondOrder][i]+1;
    279                         CreateFragment(Fragment, KeySet.OrderSet[BondOrder][i]);
    280                         for (int l=0;l<Fragment.ColumnCounter;l++)
    281                                 output << scientific << "\t" << Fragment.Matrix[ KeySet.OrderSet[BondOrder][i] ][ Fragment.RowCounter[ KeySet.OrderSet[BondOrder][i] ] ][l];
    282                         output << endl;
    283                 }
    284         }
    285         output.close();
    286         return true;
     388  stringstream filename;
     389  ofstream output;
     390
     391  filename << prefix << ".dat";
     392  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     393  cout << msg << endl;
     394  output << "# " << msg << ", created on " << datum << endl;
     395  output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
     396  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     397    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
     398      output << BondOrder+1 << "\t" << KeySet.OrderSet[BondOrder][i]+1;
     399      CreateFragment(Fragment, KeySet.OrderSet[BondOrder][i]);
     400      for (int l=0;l<Fragment.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];l++)
     401        output << scientific << "\t" << Fragment.Matrix[ KeySet.OrderSet[BondOrder][i] ][ Fragment.RowCounter[ KeySet.OrderSet[BondOrder][i] ] ][l];
     402      output << endl;
     403    }
     404  }
     405  output.close();
     406  return true;
    287407};
    288408
    … …  
    291411 * \param &KeySet KeySetsContainer with associations of each fragment to a bond order
    292412 * \param BondOrder current bond order
    293  */
     413 */ 
    294414void CreateMaxFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySet, int BondOrder)
    295415{
    296         for(int j=Fragments.RowCounter[ Fragments.MatrixCounter ];j--;) {
    297                 for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
    298                         if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
    299                                 for (int k=Fragments.ColumnCounter;k--;)
    300                                         Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    301                         }
    302                 }
    303         }
     416  for(int j=Fragments.RowCounter[ Fragments.MatrixCounter ];j--;) {
     417    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
     418      if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
     419        for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
     420          Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     421      }
     422    }     
     423  }
    304424};
    305425
    … …  
    308428 * \param &KeySet KeySetsContainer with associations of each fragment to a bond order
    309429 * \param BondOrder current bond order
    310  */
     430 */ 
    311431void CreateMinFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySet, int BondOrder)
    312432{
    313         for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
    314                 int i=0;
    315                 do {    // first get a minimum value unequal to 0
    316                         for (int k=Fragments.ColumnCounter;k--;)
    317                                 Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    318                         i++;
    319                 } while ((fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < MYEPSILON) && (i<KeySet.FragmentsPerOrder[BondOrder]));
    320                 for(;i<KeySet.FragmentsPerOrder[BondOrder];i++) { // then find lowest
    321                         if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) > fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
    322                                 for (int k=Fragments.ColumnCounter;k--;)
    323                                         Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    324                         }
    325                 }
    326         }
     433  for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
     434    int i=0;
     435    do {  // first get a minimum value unequal to 0
     436      for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
     437        Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     438      i++;
     439    } while ((fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < MYEPSILON) && (i<KeySet.FragmentsPerOrder[BondOrder]));
     440    for(;i<KeySet.FragmentsPerOrder[BondOrder];i++) { // then find lowest
     441      if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) > fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
     442        for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
     443          Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     444      }
     445    }     
     446  }
    327447};
    328448
    … …  
    331451bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int))
    332452{
    333         stringstream filename;
    334         ofstream output;
    335 
    336         filename << prefix << ".dat";
    337         if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
    338         cout << msg << endl;
    339         output << "# " << msg << ", created on " << datum;
    340         output << "#Order\tFrag.No.\t" << Fragment.Header << endl;
    341         // max
    342         for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    343                 Fragment.SetLastMatrix(0.,0);
    344                 CreateFragmentOrder(Fragment, KeySet, BondOrder);
    345                 output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    346                 for (int l=0;l<Fragment.ColumnCounter;l++)
    347                         output << scientific << "\t" << Fragment.Matrix[ Fragment.MatrixCounter ][ Fragment.RowCounter[ Fragment.MatrixCounter ]-1 ][l];
    348                 output << endl;
    349         }
    350         output.close();
    351         return true;
     453  stringstream filename;
     454  ofstream output;
     455
     456  filename << prefix << ".dat";
     457  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     458  cout << msg << endl;
     459  output << "# " << msg << ", created on " << datum;
     460  output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
     461  // max
     462  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     463    Fragment.SetLastMatrix(0.,0);
     464    CreateFragmentOrder(Fragment, KeySet, BondOrder);
     465    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     466    for (int l=0;l<Fragment.ColumnCounter[ Fragment.MatrixCounter ];l++)
     467      output << scientific << "\t" << Fragment.Matrix[ Fragment.MatrixCounter ][ Fragment.RowCounter[ Fragment.MatrixCounter ]-1 ][l];
     468    output << endl;
     469  }
     470  output.close();
     471  return true;
    352472};
    353473
    … …  
    358478void CreateEnergy(class MatrixContainer &Energy, int MatrixNumber)
    359479{
    360         for(int k=0;k<Energy.ColumnCounter;k++)
    361                 Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber] ] [k] =    Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber]-1 ] [k];
     480  for(int k=0;k<Energy.ColumnCounter[MatrixNumber];k++)
     481    Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber] ] [k] =  Energy.Matrix[MatrixNumber][ Energy.RowCounter[MatrixNumber]-1 ] [k];
    362482};
    363483
    … …  
    369489void CreateMinimumForce(class MatrixContainer &Force, int MatrixNumber)
    370490{
    371         for (int l=Force.ColumnCounter;l--;)
    372                 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
    373         for (int l=5;l<Force.ColumnCounter;l+=3) {
    374                 double stored = 0;
    375                 int k=0;
    376                 do {
    377                         for (int m=NDIM;m--;) {
    378                                 stored += Force.Matrix[MatrixNumber][ k ][l+m]
    379                                                         * Force.Matrix[MatrixNumber][ k ][l+m];
    380                                 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]       = Force.Matrix[MatrixNumber][ k ][l+m];
    381                         }
    382                         k++;
    383                 } while ((fabs(stored) < MYEPSILON) && (k<Force.RowCounter[MatrixNumber]));
    384                 for (;k<Force.RowCounter[MatrixNumber];k++) {
    385                         double tmp = 0;
    386                         for (int m=NDIM;m--;)
    387                                 tmp += Force.Matrix[MatrixNumber][ k ][l+m]
    388                                                         * Force.Matrix[MatrixNumber][ k ][l+m];
    389                         if ((fabs(tmp) > MYEPSILON) && (tmp < stored)) {        // current force is greater than stored
    390                                 for (int m=NDIM;m--;)
    391                                         Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]       = Force.Matrix[MatrixNumber][ k ][l+m];
    392                                 stored = tmp;
    393                         }
    394                 }
    395         }
     491  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
     492    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
     493  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
     494    double stored = 0;
     495    int k=0;
     496    do {
     497      for (int m=NDIM;m--;) {
     498        stored += Force.Matrix[MatrixNumber][ k ][l+m]
     499              * Force.Matrix[MatrixNumber][ k ][l+m];
     500        Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]  = Force.Matrix[MatrixNumber][ k ][l+m];
     501      }
     502      k++;
     503    } while ((fabs(stored) < MYEPSILON) && (k<Force.RowCounter[MatrixNumber]));
     504    for (;k<Force.RowCounter[MatrixNumber];k++) {
     505      double tmp = 0;
     506      for (int m=NDIM;m--;)
     507        tmp += Force.Matrix[MatrixNumber][ k ][l+m]
     508              * Force.Matrix[MatrixNumber][ k ][l+m];
     509      if ((fabs(tmp) > MYEPSILON) && (tmp < stored)) {  // current force is greater than stored
     510        for (int m=NDIM;m--;)
     511          Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]  = Force.Matrix[MatrixNumber][ k ][l+m];
     512        stored = tmp;
     513      }
     514    }
     515  }
    396516};
    397517
    … …  
    399519 * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
    400520 * \param Force ForceMatrix class containing matrix values
    401         * \param MatrixNumber the index for the ForceMatrix::matrix array
     521  * \param MatrixNumber the index for the ForceMatrix::matrix array
    402522 */
    403523void CreateMeanForce(class MatrixContainer &Force, int MatrixNumber)
    404524{
    405         int divisor = 0;
    406         for (int l=Force.ColumnCounter;l--;)
    407                 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
    408         for (int l=5;l<Force.ColumnCounter;l+=3) {
    409                 double tmp = 0;
    410                 for (int k=Force.RowCounter[MatrixNumber];k--;) {
    411                         double norm = 0.;
    412                         for (int m=NDIM;m--;)
    413                                 norm += Force.Matrix[MatrixNumber][ k ][l+m]
    414                                                         * Force.Matrix[MatrixNumber][ k ][l+m];
    415                         tmp += sqrt(norm);
    416                         if (fabs(norm) > MYEPSILON) divisor++;
    417                 }
    418                 tmp /= (double)divisor;
    419                 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = tmp;
    420         }
     525  int divisor = 0;
     526  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
     527    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
     528  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
     529    double tmp = 0;
     530    for (int k=Force.RowCounter[MatrixNumber];k--;) {
     531      double norm = 0.;
     532      for (int m=NDIM;m--;)
     533        norm += Force.Matrix[MatrixNumber][ k ][l+m]
     534              * Force.Matrix[MatrixNumber][ k ][l+m];
     535      tmp += sqrt(norm);
     536      if (fabs(norm) > MYEPSILON) divisor++;
     537    }
     538    tmp /= (double)divisor;
     539    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = tmp;
     540  }
    421541};
    422542
    … …  
    428548void CreateMaximumForce(class MatrixContainer &Force, int MatrixNumber)
    429549{
    430         for (int l=5;l<Force.ColumnCounter;l+=3) {
    431                 double stored = 0;
    432                 for (int k=Force.RowCounter[MatrixNumber];k--;) {
    433                         double tmp = 0;
    434                         for (int m=NDIM;m--;)
    435                                 tmp += Force.Matrix[MatrixNumber][ k ][l+m]
    436                                                         * Force.Matrix[MatrixNumber][ k ][l+m];
    437                         if (tmp > stored) {     // current force is greater than stored
    438                                 for (int m=NDIM;m--;)
    439                                         Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]       = Force.Matrix[MatrixNumber][ k ][l+m];
    440                                 stored = tmp;
    441                         }
    442                 }
    443         }
     550  for (int l=5;l<Force.ColumnCounter[MatrixNumber];l+=3) {
     551    double stored = 0;
     552    for (int k=Force.RowCounter[MatrixNumber];k--;) {
     553      double tmp = 0;
     554      for (int m=NDIM;m--;)
     555        tmp += Force.Matrix[MatrixNumber][ k ][l+m]
     556              * Force.Matrix[MatrixNumber][ k ][l+m];
     557      if (tmp > stored) {  // current force is greater than stored
     558        for (int m=NDIM;m--;)
     559          Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m]  = Force.Matrix[MatrixNumber][ k ][l+m];
     560        stored = tmp;
     561      }
     562    }
     563  }
    444564};
    445565
    … …  
    450570void CreateSameForce(class MatrixContainer &Force, int MatrixNumber)
    451571{
    452         // does nothing
     572  // does nothing
    453573};
    454574
    … …  
    460580void CreateVectorSumForce(class MatrixContainer &Force, int MatrixNumber)
    461581{
    462         for (int l=Force.ColumnCounter;l--;)
    463                 Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
    464         for (int l=0;l<Force.ColumnCounter;l++) {
    465                 for (int k=Force.RowCounter[MatrixNumber];k--;)
    466                         Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] += Force.Matrix[MatrixNumber][k][l];
    467         }
     582  for (int l=Force.ColumnCounter[MatrixNumber];l--;)
     583    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
     584  for (int l=0;l<Force.ColumnCounter[MatrixNumber];l++) {
     585    for (int k=Force.RowCounter[MatrixNumber];k--;)
     586      Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] += Force.Matrix[MatrixNumber][k][l];
     587  }
    468588};
    469589
    … …  
    472592 * \param *key position of key
    473593 * \param *logscale axis for logscale
    474  * \param *extraline extra set lines if desired
     594 * \param *extraline extra set lines if desired 
    475595 * \param mxtics small tics at ...
    476596 * \param xtics large tics at ...
    477  * \param *xlabel label for x axis
     597 * \param *xlabel label for x axis 
    478598 * \param *ylabel label for y axis
    479  */
     599 */ 
    480600void CreatePlotHeader(ofstream &output, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel)
    481601{
    482         //output << "#!/home/heber/build/pyxplot/pyxplot" << endl << endl;
    483         output << "reset" << endl;
    484         output << "set keycolumns "<< keycolumns << endl;
    485         output << "set key " << key << endl;
    486         output << "set mxtics "<< mxtics << endl;
    487         output << "set xtics "<< xtics << endl;
    488         if (logscale != NULL)
    489                 output << "set logscale " << logscale << endl;
    490         if (extraline != NULL)
    491                 output << extraline << endl;
    492         output << "set xlabel '" << xlabel << "'" << endl;
    493         output << "set ylabel '" << ylabel << "'" << endl;
    494         output << "set terminal eps color" << endl;
    495         output << "set output '"<< prefix << ".eps'" << endl;
     602  //output << "#!/home/heber/build/pyxplot/pyxplot" << endl << endl;
     603  output << "reset" << endl;
     604  output << "set keycolumns "<< keycolumns << endl;
     605  output << "set key " << key << endl;
     606  output << "set mxtics "<< mxtics << endl;
     607  output << "set xtics "<< xtics << endl;
     608  if (logscale != NULL)
     609    output << "set logscale " << logscale << endl;
     610  if (extraline != NULL)
     611    output << extraline << endl;
     612  output << "set xlabel '" << xlabel << "'" << endl;
     613  output << "set ylabel '" << ylabel << "'" << endl;
     614  output << "set terminal eps color" << endl;
     615  output << "set output '"<< prefix << ".eps'" << endl;
    496616};
    497617
    … …  
    506626 * \param mxtics small tics at ...
    507627 * \param xtics large tics at ...
    508  * \param xlabel label for x axis
     628 * \param xlabel label for x axis 
    509629 * \param xlabel label for x axis
    510630 * \param *xrange xrange
    … …  
    514634 * \param (*CreatePlotLines) function reference that writes a single plot line
    515635 * \return true if file was written successfully
    516  */
     636 */   
    517637bool CreatePlotOrder(class MatrixContainer &Matrix, const class KeySetsContainer &KeySet, const char *dir, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel, const char *xrange, const char *yrange, const char *xargument, const char *uses, void (*CreatePlotLines)(ofstream &, class MatrixContainer &, const char *, const char *, const char *))
    518638{
    519         stringstream filename;
    520         ofstream output;
    521 
    522         filename << prefix << ".pyx";
    523         if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
    524         CreatePlotHeader(output, prefix, keycolumns, key, logscale, extraline, mxtics, xtics, xlabel, ylabel);
    525         output << "plot " << xrange << " " << yrange << " \\" << endl;
    526         CreatePlotLines(output, Matrix, prefix, xargument, uses);
    527         output.close();
    528         return true;
     639  stringstream filename;
     640  ofstream output;
     641
     642  filename << prefix << ".pyx";
     643  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     644  CreatePlotHeader(output, prefix, keycolumns, key, logscale, extraline, mxtics, xtics, xlabel, ylabel);
     645  output << "plot " << xrange << " " << yrange << " \\" << endl;
     646  CreatePlotLines(output, Matrix, prefix, xargument, uses);
     647  output.close(); 
     648  return true;
    529649};
    530650
    … …  
    538658void AbsEnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
    539659{
    540         stringstream line(Energy.Header);
    541         string token;
    542 
    543         getline(line, token, '\t');
    544         for (int i=2; i<= Energy.ColumnCounter;i++) {
    545                 getline(line, token, '\t');
    546                 while (token[0] == ' ') // remove leading white spaces
    547                         token.erase(0,1);
    548                 output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(abs($" << i+2 << ")) " << uses;
    549                 if (i != (Energy.ColumnCounter))
    550                         output << ", \\";
    551                 output << endl;
    552         }
     660  stringstream line(Energy.Header[ Energy.MatrixCounter ]);
     661  string token;
     662
     663  getline(line, token, '\t');
     664  for (int i=2; i<= Energy.ColumnCounter[Energy.MatrixCounter];i++) {
     665    getline(line, token, '\t');
     666    while (token[0] == ' ') // remove leading white spaces
     667      token.erase(0,1);
     668    output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(abs($" << i+2 << ")) " << uses;
     669    if (i != (Energy.ColumnCounter[Energy.MatrixCounter]))
     670      output << ", \\";
     671    output << endl;
     672  }
    553673};
    554674
    … …  
    562682void EnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
    563683{
    564         stringstream line(Energy.Header);
    565         string token;
    566 
    567         getline(line, token, '\t');
    568         for (int i=1; i<= Energy.ColumnCounter;i++) {
    569                 getline(line, token, '\t');
    570                 while (token[0] == ' ') // remove leading white spaces
    571                         token.erase(0,1);
    572                 output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":" << i+2 << " " << uses;
    573                 if (i != (Energy.ColumnCounter))
    574                         output << ", \\";
    575                 output << endl;
    576         }
     684  stringstream line(Energy.Header[Energy.MatrixCounter]);
     685  string token;
     686
     687  getline(line, token, '\t');
     688  for (int i=1; i<= Energy.ColumnCounter[Energy.MatrixCounter];i++) {
     689    getline(line, token, '\t');
     690    while (token[0] == ' ') // remove leading white spaces
     691      token.erase(0,1);
     692    output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":" << i+2 << " " << uses;
     693    if (i != (Energy.ColumnCounter[Energy.MatrixCounter]))
     694      output << ", \\";
     695    output << endl;
     696  }
    577697};
    578698
    … …  
    586706void ForceMagnitudePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
    587707{
    588         stringstream line(Force.Header);
    589         string token;
    590 
    591         getline(line, token, '\t');
    592         getline(line, token, '\t');
    593         getline(line, token, '\t');
    594         getline(line, token, '\t');
    595         getline(line, token, '\t');
    596         for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
    597                 getline(line, token, '\t');
    598                 while (token[0] == ' ') // remove leading white spaces
    599                         token.erase(0,1);
    600                 token.erase(token.length(), 1); // kill residual index char (the '0')
    601                 output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses;
    602                 if (i != (Force.ColumnCounter-1))
    603                         output << ", \\";
    604                 output << endl;
    605                 getline(line, token, '\t');
    606                 getline(line, token, '\t');
    607         }
     708  stringstream line(Force.Header[Force.MatrixCounter]);
     709  string token;
     710
     711  getline(line, token, '\t');
     712  getline(line, token, '\t');
     713  getline(line, token, '\t');
     714  getline(line, token, '\t');
     715  getline(line, token, '\t');
     716  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     717    getline(line, token, '\t');
     718    while (token[0] == ' ') // remove leading white spaces
     719      token.erase(0,1);
     720    token.erase(token.length(), 1);  // kill residual index char (the '0')
     721    output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses;
     722    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
     723      output << ", \\";
     724    output << endl;
     725    getline(line, token, '\t');
     726    getline(line, token, '\t');
     727  }
    608728};
    609729
    … …  
    617737void AbsFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
    618738{
    619         stringstream line(Force.Header);
    620         string token;
    621 
    622         getline(line, token, '\t');
    623         getline(line, token, '\t');
    624         getline(line, token, '\t');
    625         getline(line, token, '\t');
    626         getline(line, token, '\t');
    627         for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
    628                 getline(line, token, '\t');
    629                 while (token[0] == ' ') // remove leading white spaces
    630                         token.erase(0,1);
    631                 token.erase(token.length(), 1); // kill residual index char (the '0')
    632                 output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
    633                 if (i != (Force.ColumnCounter-1))
    634                         output << ", \\";
    635                 output << endl;
    636                 getline(line, token, '\t');
    637                 getline(line, token, '\t');
    638         }
     739  stringstream line(Force.Header[Force.MatrixCounter]);
     740  string token;
     741
     742  getline(line, token, '\t');
     743  getline(line, token, '\t');
     744  getline(line, token, '\t');
     745  getline(line, token, '\t');
     746  getline(line, token, '\t');
     747  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     748    getline(line, token, '\t');
     749    while (token[0] == ' ') // remove leading white spaces
     750      token.erase(0,1);
     751    token.erase(token.length(), 1);  // kill residual index char (the '0')
     752    output << "'" << prefix << ".dat' title '" << token << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
     753    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
     754      output << ", \\";
     755    output << endl;
     756    getline(line, token, '\t');
     757    getline(line, token, '\t');
     758  }
    639759};
    640760
    … …  
    648768void BoxesForcePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
    649769{
    650         stringstream line(Force.Header);
    651         char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
    652         string token;
    653 
    654         getline(line, token, '\t');
    655         getline(line, token, '\t');
    656         getline(line, token, '\t');
    657         getline(line, token, '\t');
    658         getline(line, token, '\t');
    659         for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
    660                 getline(line, token, '\t');
    661                 while (token[0] == ' ') // remove leading white spaces
    662                         token.erase(0,1);
    663                 token.erase(token.length(), 1); // kill residual index char (the '0')
    664                 output << "'" << prefix << ".dat' title '" << token << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses << " " << fillcolor[(i-7)/3];
    665                 if (i != (Force.ColumnCounter-1))
    666                         output << ", \\";
    667                 output << endl;
    668                 getline(line, token, '\t');
    669                 getline(line, token, '\t');
    670         }
     770  stringstream line(Force.Header[Force.MatrixCounter]);
     771  char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
     772  string token;
     773
     774  getline(line, token, '\t');
     775  getline(line, token, '\t');
     776  getline(line, token, '\t');
     777  getline(line, token, '\t');
     778  getline(line, token, '\t');
     779  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     780    getline(line, token, '\t');
     781    while (token[0] == ' ') // remove leading white spaces
     782      token.erase(0,1);
     783    token.erase(token.length(), 1);  // kill residual index char (the '0')
     784    output << "'" << prefix << ".dat' title '" << token << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):(sqrt($" << i+1 << "*$" << i+1 << "+$" << i+2 << "*$" << i+2 << "+$" << i+3 << "*$" << i+3 << ")) " << uses << " " << fillcolor[(i-7)/3];
     785    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
     786      output << ", \\";
     787    output << endl;
     788    getline(line, token, '\t');
     789    getline(line, token, '\t');
     790  }
    671791};
    672792
    … …  
    680800void BoxesFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
    681801{
    682         stringstream line(Force.Header);
    683         char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
    684         string token;
    685 
    686         getline(line, token, '\t');
    687         getline(line, token, '\t');
    688         getline(line, token, '\t');
    689         getline(line, token, '\t');
    690         getline(line, token, '\t');
    691         for (int i=7; i< Force.ColumnCounter;i+=NDIM) {
    692                 getline(line, token, '\t');
    693                 while (token[0] == ' ') // remove leading white spaces
    694                         token.erase(0,1);
    695                 token.erase(token.length(), 1); // kill residual index char (the '0')
    696                 output << "'" << prefix << ".dat' title '" << token << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):" << i+1 << " " << uses << " " << fillcolor[(i-7)/3];
    697                 if (i != (Force.ColumnCounter-1))
    698                         output << ", \\";
    699                 output << endl;
    700                 getline(line, token, '\t');
    701                 getline(line, token, '\t');
    702         }
    703 };
     802  stringstream line(Force.Header[Force.MatrixCounter]);
     803  char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
     804  string token;
     805
     806  getline(line, token, '\t');
     807  getline(line, token, '\t');
     808  getline(line, token, '\t');
     809  getline(line, token, '\t');
     810  getline(line, token, '\t');
     811  for (int i=7; i< Force.ColumnCounter[Force.MatrixCounter];i+=NDIM) {
     812    getline(line, token, '\t');
     813    while (token[0] == ' ') // remove leading white spaces
     814      token.erase(0,1);
     815    token.erase(token.length(), 1);  // kill residual index char (the '0')
     816    output << "'" << prefix << ".dat' title '" << token << "' using ($" << xargument << "+" << fixed << setprecision(1) << (double)((i-7)/3)*0.2 << "):" << i+1 << " " << uses << " " << fillcolor[(i-7)/3];
     817    if (i != (Force.ColumnCounter[Force.MatrixCounter]-1))
     818      output << ", \\";
     819    output << endl;
     820    getline(line, token, '\t');
     821    getline(line, token, '\t');
     822  }
     823};
Note: See TracChangeset for help on using the changeset viewer.