Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/parser.cpp

    rce5ac3 r6ac7ee  
    22 *
    33 * Declarations of various class functions for the parsing of value files.
    4  *   
     4 *
    55 */
    66
    77// ======================================= INCLUDES ==========================================
    88
    9 #include "helpers.hpp" 
     9#include "helpers.hpp"
    1010#include "parser.hpp"
    1111
    … …  
    2424bool FilePresent(const char *filename, bool test)
    2525{
    26   ifstream input;
    27  
    28   input.open(filename, ios::in);
    29   if (input == NULL) {
    30     if (!test)
    31       cout << endl << "Unable to open " << filename << ", is the directory correct?" << endl;
    32     return false;
    33   }
    34   input.close();
    35   return true;
     26        ifstream input;
     27
     28        input.open(filename, ios::in);
     29        if (input == NULL) {
     30                if (!test)
     31                        cout << endl << "Unable to open " << filename << ", is the directory correct?" << endl;
     32                return false;
     33        }
     34        input.close();
     35        return true;
    3636};
    3737
    … …  
    4343bool TestParams(int argc, char **argv)
    4444{
    45   ifstream input;
    46   stringstream line;
    47 
    48   line << argv[1] << FRAGMENTPREFIX << KEYSETFILE;
    49   return FilePresent(line.str().c_str(), false);
     45        ifstream input;
     46        stringstream line;
     47
     48        line << argv[1] << FRAGMENTPREFIX << KEYSETFILE;
     49        return FilePresent(line.str().c_str(), false);
    5050};
    5151
    … …  
    5555 */
    5656MatrixContainer::MatrixContainer() {
    57   Indices = NULL;
    58   Header = (char **) Malloc(sizeof(char)*1, "MatrixContainer::MatrixContainer: **Header");
    59   Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
    60   RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter");
    61   ColumnCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *ColumnCounter");
    62   Header[0] = NULL;
    63   Matrix[0] = NULL;
    64   RowCounter[0] = -1;
    65   MatrixCounter = 0;
    66   ColumnCounter[0] = -1;
     57        Indices = NULL;
     58        Header = (char *) Malloc(sizeof(char)*1023, "MatrixContainer::MatrixContainer: *Header");
     59        Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
     60        RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter");
     61        Matrix[0] = NULL;
     62        RowCounter[0] = -1;
     63        MatrixCounter = 0;
     64        ColumnCounter = 0;
    6765};
    6866
    … …  
    7068 */
    7169MatrixContainer::~MatrixContainer() {
    72   if (Matrix != NULL) {
    73     for(int i=MatrixCounter;i--;) {
    74       if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
    75           for(int j=RowCounter[i]+1;j--;)
    76             Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
    77         Free((void **)&Matrix[i], "MatrixContainer::~MatrixContainer: **Matrix[]");
    78       }
    79     }
    80     if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
    81       for(int j=RowCounter[MatrixCounter]+1;j--;)
    82         Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
    83     if (MatrixCounter != 0)
    84       Free((void **)&Matrix[MatrixCounter], "MatrixContainer::~MatrixContainer: **Matrix[MatrixCounter]");
    85     Free((void **)&Matrix, "MatrixContainer::~MatrixContainer: ***Matrix");
    86   }
    87   if (Indices != NULL)
    88     for(int i=MatrixCounter+1;i--;) {
    89       Free((void **)&Indices[i], "MatrixContainer::~MatrixContainer: *Indices[]");
    90     }
    91   Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices");
    92  
    93   if (Header != NULL)
    94     for(int i=MatrixCounter+1;i--;)
    95       Free((void **)&Header[i], "MatrixContainer::~MatrixContainer: *Header[]");
    96   Free((void **)&Header, "MatrixContainer::~MatrixContainer: **Header");
    97   Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
    98   Free((void **)&ColumnCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
    99 };
    100 
    101 /** Either copies index matrix from another MatrixContainer or initialises with trivial mapping if NULL.
    102  * This either copies the index matrix or just maps 1 to 1, 2 to 2 and so on for all fragments.
    103  * \param *Matrix pointer to other MatrixContainer
    104  * \return true - copy/initialisation sucessful, false - dimension false for copying
    105  */
    106 bool MatrixContainer::InitialiseIndices(class MatrixContainer *Matrix)
    107 {
    108   cout << "Initialising indices";
    109   if (Matrix == NULL) {
    110     cout << " with trivial mapping." << endl;
    111     Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "MatrixContainer::InitialiseIndices: **Indices");
    112     for(int i=MatrixCounter+1;i--;) {
    113       Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
    114       for(int j=RowCounter[i];j--;)
    115         Indices[i][j] = j;
    116     }
    117   } else {
    118     cout << " from other MatrixContainer." << endl;
    119     if (MatrixCounter != Matrix->MatrixCounter)
    120       return false;
    121     Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "MatrixContainer::InitialiseIndices: **Indices");
    122     for(int i=MatrixCounter+1;i--;) {
    123       if (RowCounter[i] != Matrix->RowCounter[i])
    124         return false;
    125       Indices[i] = (int *) Malloc(sizeof(int)*Matrix->RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
    126       for(int j=Matrix->RowCounter[i];j--;) {
    127         Indices[i][j] = Matrix->Indices[i][j];
    128         //cout << Indices[i][j] << "\t";
    129       }
    130       //cout << endl;
    131     }
    132   }
    133   return true;
    134 };
     70        if (Matrix != NULL) {
     71                for(int i=MatrixCounter;i--;) {
     72                        if (RowCounter != NULL) {
     73                                        for(int j=RowCounter[i]+1;j--;)
     74                                                Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
     75                                Free((void **)&Matrix[i], "MatrixContainer::~MatrixContainer: **Matrix[]");
     76                        }
     77                }
     78                if ((RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
     79                        for(int j=RowCounter[MatrixCounter]+1;j--;)
     80                                Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
     81                if (MatrixCounter != 0)
     82                        Free((void **)&Matrix[MatrixCounter], "MatrixContainer::~MatrixContainer: **Matrix[MatrixCounter]");
     83                Free((void **)&Matrix, "MatrixContainer::~MatrixContainer: ***Matrix");
     84        }
     85        if (Indices != NULL)
     86                for(int i=MatrixCounter+1;i--;) {
     87                        Free((void **)&Indices[i], "MatrixContainer::~MatrixContainer: *Indices[]");
     88                }
     89        Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices");
     90
     91        Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
     92        Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
     93};
     94
    13595
    13696/** Parsing a number of matrices.
    137  *    -# open the matrix file
    138  *    -# skip some lines (\a skiplines)
    139  *    -# scan header lines for number of columns
    140  *    -# scan lines for number of rows
    141  *    -# allocate matrix
    142  *    -# loop over found column and row counts and parse in each entry
     97 *              -# open the matrix file
     98 *              -# skip some lines (\a skiplines)
     99 *              -# scan header lines for number of columns
     100 *              -# scan lines for number of rows
     101 *              -# allocate matrix
     102 *              -# loop over found column and row counts and parse in each entry
    143103 * \param *name directory with files
    144104 * \param skiplines number of inital lines to skip
    145105 * \param skiplines number of inital columns to skip
    146106 * \param MatrixNr index number in Matrix array to parse into
    147  * \return parsing successful 
     107 * \return parsing successful
    148108 */
    149109bool MatrixContainer::ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr)
    150110{
    151   ifstream input;
    152   stringstream line;
    153   string token;
    154   char filename[1023];
    155 
    156   input.open(name, ios::in);
    157   //cout << "Opening " << name << " ... "  << input << endl;
    158   if (input == NULL) {
    159     cerr << endl << "Unable to open " << name << ", is the directory correct?" << endl;
    160     return false;
    161   }
    162 
    163   // parse header
    164   Header[MatrixNr] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseMatrix: *Header[]");
    165   for (int m=skiplines+1;m--;)
    166     input.getline(Header[MatrixNr], 1023);
    167  
    168   // scan header for number of columns
    169   line.str(Header[MatrixNr]);
    170   for(int k=skipcolumns;k--;)
    171     line >> Header[MatrixNr];
    172   //cout << line.str() << endl;
    173   ColumnCounter[MatrixNr]=0;
    174   while ( getline(line,token, '\t') ) {
    175     if (token.length() > 0)
    176       ColumnCounter[MatrixNr]++;
    177   }
    178   //cout << line.str() << endl;
    179   //cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
    180   if (ColumnCounter[MatrixNr] == 0)
    181     cerr << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
    182  
    183   // scan rest for number of rows/lines
    184   RowCounter[MatrixNr]=-1;    // counts one line too much
    185   while (!input.eof()) {
    186     input.getline(filename, 1023);
    187     //cout << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
    188     RowCounter[MatrixNr]++; // then line was not last MeanForce
    189     if (strncmp(filename,"MeanForce", 9) == 0) {
    190       break;
    191     }
    192   }
    193   //cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
    194   if (RowCounter[MatrixNr] == 0)
    195     cerr << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
    196  
    197   // allocate matrix if it's not zero dimension in one direction
    198   if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
    199     Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseMatrix: **Matrix[]");
    200  
    201     // parse in each entry for this matrix
    202     input.clear();
    203     input.seekg(ios::beg);
    204     for (int m=skiplines+1;m--;)
    205       input.getline(Header[MatrixNr], 1023);    // skip header
    206     line.str(Header[MatrixNr]);
    207     for(int k=skipcolumns;k--;)  // skip columns in header too
    208       line >> filename;
    209     strncpy(Header[MatrixNr], line.str().c_str(), 1023); 
    210     for(int j=0;j<RowCounter[MatrixNr];j++) {
    211       Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]");
    212       input.getline(filename, 1023);
    213       stringstream lines(filename);
    214       //cout << "Matrix at level " << j << ":";// << filename << endl;
    215       for(int k=skipcolumns;k--;)
    216         lines >> filename;
    217       for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
    218         lines >> Matrix[MatrixNr][j][k];
    219         //cout << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
    220       }
    221       //cout << endl;
    222       Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]");
    223       for(int j=ColumnCounter[MatrixNr];j--;)
    224         Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
    225     }
    226   } else {
    227     cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl;
    228   }
    229   input.close();
    230   return true;
     111        ifstream input;
     112        stringstream line;
     113        string token;
     114        char filename[1023];
     115
     116        input.open(name, ios::in);
     117        //cout << "Opening " << name << " ... " << input << endl;
     118        if (input == NULL) {
     119                cerr << endl << "Unable to open " << name << ", is the directory correct?" << endl;
     120                return false;
     121        }
     122
     123        // skip some initial lines
     124        for (int m=skiplines+1;m--;)
     125                input.getline(Header, 1023);
     126
     127        // scan header for number of columns
     128        line.str(Header);
     129        for(int k=skipcolumns;k--;)
     130                line >> Header;
     131        //cout << line.str() << endl;
     132        ColumnCounter=0;
     133        while ( getline(line,token, '\t') ) {
     134                if (token.length() > 0)
     135                        ColumnCounter++;
     136        }
     137        //cout << line.str() << endl;
     138        //cout << "ColumnCounter: " << ColumnCounter << "." << endl;
     139        if (ColumnCounter == 0)
     140                cerr << "ColumnCounter: " << ColumnCounter << " from file " << name << ", this is probably an error!" << endl;
     141
     142        // scan rest for number of rows/lines
     143        RowCounter[MatrixNr]=-1;                // counts one line too much
     144        while (!input.eof()) {
     145                input.getline(filename, 1023);
     146                //cout << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
     147                RowCounter[MatrixNr]++; // then line was not last MeanForce
     148                if (strncmp(filename,"MeanForce", 9) == 0) {
     149                        break;
     150                }
     151        }
     152        //cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
     153        if (RowCounter[MatrixNr] == 0)
     154                cerr << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
     155
     156        // allocate matrix if it's not zero dimension in one direction
     157        if ((ColumnCounter > 0) && (RowCounter[MatrixNr] > -1)) {
     158                Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     159
     160                // parse in each entry for this matrix
     161                input.clear();
     162                input.seekg(ios::beg);
     163                for (int m=skiplines+1;m--;)
     164                        input.getline(Header, 1023);            // skip header
     165                line.str(Header);
     166                for(int k=skipcolumns;k--;)     // skip columns in header too
     167                        line >> filename;
     168                strncpy(Header, line.str().c_str(), 1023);
     169                for(int j=0;j<RowCounter[MatrixNr];j++) {
     170                        Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     171                        input.getline(filename, 1023);
     172                        stringstream lines(filename);
     173                        //cout << "Matrix at level " << j << ":";// << filename << endl;
     174                        for(int k=skipcolumns;k--;)
     175                                lines >> filename;
     176                        for(int k=0;(k<ColumnCounter) && (!lines.eof());k++) {
     177                                lines >> Matrix[MatrixNr][j][k];
     178                                //cout << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
     179                        }
     180                        //cout << endl;
     181                        Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
     182                        for(int j=ColumnCounter;j--;)
     183                                Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
     184                }
     185        } else {
     186                cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl;
     187        }
     188        input.close();
     189        return true;
    231190};
    232191
    233192/** Parsing a number of matrices.
    234193 * -# First, count the number of matrices by counting lines in KEYSETFILE
    235  * -# Then, 
    236  *    -# construct the fragment number
    237  *    -# open the matrix file
    238  *    -# skip some lines (\a skiplines)
    239  *    -# scan header lines for number of columns
    240  *    -# scan lines for number of rows
    241  *    -# allocate matrix
    242  *    -# loop over found column and row counts and parse in each entry
    243  * -# Finally, allocate one additional matrix (\a MatrixCounter) containing combined or temporary values 
     194 * -# Then,
     195 *              -# construct the fragment number
     196 *              -# open the matrix file
     197 *              -# skip some lines (\a skiplines)
     198 *              -# scan header lines for number of columns
     199 *              -# scan lines for number of rows
     200 *              -# allocate matrix
     201 *              -# loop over found column and row counts and parse in each entry
     202 * -# Finally, allocate one additional matrix (\a MatrixCounter) containing combined or temporary values
    244203 * \param *name directory with files
    245204 * \param *prefix prefix of each matrix file
    … …  
    247206 * \param skiplines number of inital lines to skip
    248207 * \param skiplines number of inital columns to skip
    249  * \return parsing successful 
     208 * \return parsing successful
    250209 */
    251210bool MatrixContainer::ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns)
    252211{
    253   char filename[1023];
    254   ifstream input;
    255   char *FragmentNumber = NULL;
    256   stringstream file;
    257   string token;
    258  
    259   // count the number of matrices
    260   MatrixCounter = -1; // we count one too much
    261   file << name << FRAGMENTPREFIX << KEYSETFILE;
    262   input.open(file.str().c_str(), ios::in);
    263   if (input == NULL) {
    264     cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
    265     return false;
    266   }
    267   while (!input.eof()) {
    268     input.getline(filename, 1023);
    269     stringstream zeile(filename);
    270     MatrixCounter++;
    271   }
    272   input.close(); 
    273   cout << "Determined " << MatrixCounter << " fragments." << endl;
    274  
    275   cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
    276   Header = (char **) ReAlloc(Header, sizeof(char *)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: **Header"); // one more each for the total molecule
    277   Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
    278   RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
    279   ColumnCounter = (int *) ReAlloc(ColumnCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *ColumnCounter");
    280   for(int i=MatrixCounter+1;i--;) {
    281     Matrix[i] = NULL;
    282     Header[i] = NULL;
    283     RowCounter[i] = -1;
    284     ColumnCounter[i] = -1;
    285   }
    286   for(int i=0; i < MatrixCounter;i++) {
    287     // open matrix file
    288     FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    289     file.str(" ");
    290     file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
    291     if (!ParseMatrix(file.str().c_str(), skiplines, skipcolumns, i))
    292       return false;
    293     Free((void **)&FragmentNumber, "MatrixContainer::ParseFragmentMatrix: *FragmentNumber");
    294   }
    295   return true;
     212        char filename[1023];
     213        ifstream input;
     214        char *FragmentNumber = NULL;
     215        stringstream file;
     216        string token;
     217
     218        // count the number of matrices
     219        MatrixCounter = -1; // we count one too much
     220        file << name << FRAGMENTPREFIX << KEYSETFILE;
     221        input.open(file.str().c_str(), ios::in);
     222        if (input == NULL) {
     223                cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
     224                return false;
     225        }
     226        while (!input.eof()) {
     227                input.getline(filename, 1023);
     228                stringstream zeile(filename);
     229                MatrixCounter++;
     230        }
     231        input.close();
     232        cout << "Determined " << MatrixCounter << " fragments." << endl;
     233
     234        cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
     235        Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
     236        RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
     237        for(int i=MatrixCounter+1;i--;) {
     238                Matrix[i] = NULL;
     239                RowCounter[i] = -1;
     240        }
     241        for(int i=0; i < MatrixCounter;i++) {
     242                // open matrix file
     243                FragmentNumber = FixedDigitNumber(MatrixCounter, i);
     244                file.str(" ");
     245                file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
     246                if (!ParseMatrix(file.str().c_str(), skiplines, skipcolumns, i))
     247                        return false;
     248                Free((void **)&FragmentNumber, "MatrixContainer::ParseFragmentMatrix: *FragmentNumber");
     249        }
     250        return true;
    296251};
    297252
    298253/** Allocates and resets the memory for a number \a MCounter of matrices.
    299  * \param **GivenHeader Header line for each matrix
     254 * \param *GivenHeader Header line
    300255 * \param MCounter number of matrices
    301256 * \param *RCounter number of rows for each matrix
    302  * \param *CCounter number of columns for each matrix
    303  * \return Allocation successful
    304  */
    305 bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
    306 {
    307   MatrixCounter = MCounter;
    308   Header = (char **) Malloc(sizeof(char *)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *Header");
    309   Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: ***Matrix"); // one more each for the total molecule
    310   RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *RowCounter");
    311   ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *ColumnCounter");
    312   for(int i=MatrixCounter+1;i--;) {
    313     Header[i] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::AllocateMatrix: *Header[i]");
    314     strncpy(Header[i], GivenHeader[i], 1023);
    315     RowCounter[i] = RCounter[i];
    316     ColumnCounter[i] = CCounter[i];
    317     Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::AllocateMatrix: **Matrix[]"); 
    318     for(int j=RowCounter[i]+1;j--;) {
    319       Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]");
    320       for(int k=ColumnCounter[i];k--;)
    321         Matrix[i][j][k] = 0.;
    322     }
    323   }
    324   return true;
     257 * \param CCounter number of columns for all matrices
     258 * \return Allocation successful
     259 */
     260bool MatrixContainer::AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int CCounter)
     261{
     262        Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader");
     263        strncpy(Header, GivenHeader, 1023);
     264
     265        MatrixCounter = MCounter;
     266        ColumnCounter = CCounter;
     267        Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
     268        RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
     269        for(int i=MatrixCounter+1;i--;) {
     270                RowCounter[i] = RCounter[i];
     271                Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     272                for(int j=RowCounter[i]+1;j--;) {
     273                        Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     274                        for(int k=ColumnCounter;k--;)
     275                                Matrix[i][j][k] = 0.;
     276                }
     277        }
     278        return true;
    325279};
    326280
    … …  
    330284bool MatrixContainer::ResetMatrix()
    331285{
    332   for(int i=MatrixCounter+1;i--;)
    333     for(int j=RowCounter[i]+1;j--;)
    334       for(int k=ColumnCounter[i];k--;)
    335         Matrix[i][j][k] = 0.;
    336   return true;
     286        for(int i=MatrixCounter+1;i--;)
     287                for(int j=RowCounter[i]+1;j--;)
     288                        for(int k=ColumnCounter;k--;)
     289                                Matrix[i][j][k] = 0.;
     290        return true;
    337291};
    338292
    … …  
    342296double MatrixContainer::FindMaxValue()
    343297{
    344   double max = Matrix[0][0][0];
    345   for(int i=MatrixCounter+1;i--;)
    346     for(int j=RowCounter[i]+1;j--;)
    347       for(int k=ColumnCounter[i];k--;)
    348         if (fabs(Matrix[i][j][k]) > max)
    349           max = fabs(Matrix[i][j][k]);
    350   if (fabs(max) < MYEPSILON)
    351     max += MYEPSILON;
     298        double max = Matrix[0][0][0];
     299        for(int i=MatrixCounter+1;i--;)
     300                for(int j=RowCounter[i]+1;j--;)
     301                        for(int k=ColumnCounter;k--;)
     302                                if (fabs(Matrix[i][j][k]) > max)
     303                                        max = fabs(Matrix[i][j][k]);
     304        if (fabs(max) < MYEPSILON)
     305                max += MYEPSILON;
    352306 return max;
    353307};
    … …  
    358312double MatrixContainer::FindMinValue()
    359313{
    360   double min = Matrix[0][0][0];
    361   for(int i=MatrixCounter+1;i--;)
    362     for(int j=RowCounter[i]+1;j--;)
    363       for(int k=ColumnCounter[i];k--;)
    364         if (fabs(Matrix[i][j][k]) < min)
    365           min = fabs(Matrix[i][j][k]);
    366   if (fabs(min) < MYEPSILON)
    367     min += MYEPSILON;
    368   return min;
     314        double min = Matrix[0][0][0];
     315        for(int i=MatrixCounter+1;i--;)
     316                for(int j=RowCounter[i]+1;j--;)
     317                        for(int k=ColumnCounter;k--;)
     318                                if (fabs(Matrix[i][j][k]) < min)
     319                                        min = fabs(Matrix[i][j][k]);
     320        if (fabs(min) < MYEPSILON)
     321                min += MYEPSILON;
     322        return min;
    369323};
    370324
    … …  
    376330bool MatrixContainer::SetLastMatrix(double value, int skipcolumns)
    377331{
    378   for(int j=RowCounter[MatrixCounter]+1;j--;)
    379     for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
    380       Matrix[MatrixCounter][j][k] = value;
    381   return true;
     332        for(int j=RowCounter[MatrixCounter]+1;j--;)
     333                for(int k=skipcolumns;k<ColumnCounter;k++)
     334                        Matrix[MatrixCounter][j][k] = value;
     335        return true;
    382336};
    383337
    … …  
    389343bool MatrixContainer::SetLastMatrix(double **values, int skipcolumns)
    390344{
    391   for(int j=RowCounter[MatrixCounter]+1;j--;)
    392     for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
    393       Matrix[MatrixCounter][j][k] = values[j][k];
    394   return true;
    395 };
    396 
    397 /** Sums the entries with each factor and put into last element of \a ***Matrix.
     345        for(int j=RowCounter[MatrixCounter]+1;j--;)
     346                for(int k=skipcolumns;k<ColumnCounter;k++)
     347                        Matrix[MatrixCounter][j][k] = values[j][k];
     348        return true;
     349};
     350
     351/** Sums the energy with each factor and put into last element of \a ***Energies.
    398352 * Sums over "E"-terms to create the "F"-terms
    399  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     353 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
    400354 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    401355 * \param Order bond order
    … …  
    404358bool MatrixContainer::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
    405359{
    406   // go through each order
    407   for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
    408     //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
    409     // then go per order through each suborder and pick together all the terms that contain this fragment
    410     for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
    411       for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
    412         if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
    413           //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
    414           // if the fragment's indices are all in the current fragment
    415           for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
    416             int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
    417             //cout << "Current index is " << k << "/" << m << "." << endl;
    418             if (m != -1) { // if it's not an added hydrogen
    419               for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
    420                 //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
    421                 if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
    422                   m = l;
    423                   break; 
    424                 }
    425               }
    426               //cout << "Corresponding index in CurrentFragment is " << m << "." << endl;
    427               if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
    428                 cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]  << " current force index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
    429                 return false;
    430               }
    431               if (Order == SubOrder) { // equal order is always copy from Energies
    432                 for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column
    433                   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
    434               } else {
    435                 for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;)
    436                   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
    437               }
    438             }
    439             //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
    440              //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
    441           }
    442         } else {
    443           //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
    444         }
    445       }
    446     }
    447    //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
    448   }
    449  
    450   return true;
     360        // go through each order
     361        for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
     362                //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     363                // then go per order through each suborder and pick together all the terms that contain this fragment
     364                for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
     365                        for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
     366                                if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
     367                                        //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
     368                                        // if the fragment's indices are all in the current fragment
     369                                        for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
     370                                                int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
     371                                                //cout << "Current index is " << k << "/" << m << "." << endl;
     372                                                if (m != -1) { // if it's not an added hydrogen
     373                                                        for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
     374                                                                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
     375                                                                if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
     376                                                                        m = l;
     377                                                                        break;
     378                                                                }
     379                                                        }
     380                                                        //cout << "Corresponding index in CurrentFragment is " << m << "." << endl;
     381                                                        if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
     382                                                                cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]    << " current force index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
     383                                                                return false;
     384                                                        }
     385                                                        if (Order == SubOrder) { // equal order is always copy from Energies
     386                                                                for(int l=ColumnCounter;l--;) // then adds/subtract each column
     387                                                                        Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     388                                                        } else {
     389                                                                for(int l=ColumnCounter;l--;)
     390                                                                        Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     391                                                        }
     392                                                }
     393                                                //if ((ColumnCounter>1) && (RowCounter[0]-1 >= 1))
     394                                                 //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
     395                                        }
     396                                } else {
     397                                        //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     398                                }
     399                        }
     400                }
     401         //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<     Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
     402        }
     403
     404        return true;
    451405};
    452406
    … …  
    458412bool MatrixContainer::WriteTotalFragments(const char *name, const char *prefix)
    459413{
    460   ofstream output;
    461   char *FragmentNumber = NULL;
    462 
    463   cout << "Writing fragment files." << endl;
    464   for(int i=0;i<MatrixCounter;i++) {
    465     stringstream line;
    466     FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    467     line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
    468     Free((void **)&FragmentNumber, "*FragmentNumber");
    469     output.open(line.str().c_str(), ios::out);
    470     if (output == NULL) {
    471       cerr << "Unable to open output energy file " << line.str() << "!" << endl;
    472       return false;
    473     }
    474     output << Header[i] << endl;
    475     for(int j=0;j<RowCounter[i];j++) {
    476       for(int k=0;k<ColumnCounter[i];k++)
    477         output << scientific << Matrix[i][j][k] << "\t";
    478       output << endl;
    479     }
    480     output.close();
    481   }
    482   return true;
     414        ofstream output;
     415        char *FragmentNumber = NULL;
     416
     417        cout << "Writing fragment files." << endl;
     418        for(int i=0;i<MatrixCounter;i++) {
     419                stringstream line;
     420                FragmentNumber = FixedDigitNumber(MatrixCounter, i);
     421                line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
     422                Free((void **)&FragmentNumber, "*FragmentNumber");
     423                output.open(line.str().c_str(), ios::out);
     424                if (output == NULL) {
     425                        cerr << "Unable to open output energy file " << line.str() << "!" << endl;
     426                        return false;
     427                }
     428                output << Header << endl;
     429                for(int j=0;j<RowCounter[i];j++) {
     430                        for(int k=0;k<ColumnCounter;k++)
     431                                output << scientific << Matrix[i][j][k] << "\t";
     432                        output << endl;
     433                }
     434                output.close();
     435        }
     436        return true;
    483437};
    484438
    … …  
    491445bool MatrixContainer::WriteLastMatrix(const char *name, const char *prefix, const char *suffix)
    492446{
    493   ofstream output;
    494   stringstream line;
    495 
    496   cout << "Writing matrix values of " << suffix << "." << endl;
    497   line << name << prefix << suffix;
    498   output.open(line.str().c_str(), ios::out);
    499   if (output == NULL) {
    500     cerr << "Unable to open output matrix file " << line.str() << "!" << endl;
    501     return false;
    502   }
    503   output << Header[MatrixCounter] << endl;
    504   for(int j=0;j<RowCounter[MatrixCounter];j++) {
    505     for(int k=0;k<ColumnCounter[MatrixCounter];k++)
    506       output << scientific << Matrix[MatrixCounter][j][k] << "\t";
    507     output << endl;
    508   }
    509   output.close();
    510   return true;
     447        ofstream output;
     448        stringstream line;
     449
     450        cout << "Writing matrix values of " << suffix << "." << endl;
     451        line << name << prefix << suffix;
     452        output.open(line.str().c_str(), ios::out);
     453        if (output == NULL) {
     454                cerr << "Unable to open output matrix file " << line.str() << "!" << endl;
     455                return false;
     456        }
     457        output << Header << endl;
     458        for(int j=0;j<RowCounter[MatrixCounter];j++) {
     459                for(int k=0;k<ColumnCounter;k++)
     460                        output << scientific << Matrix[MatrixCounter][j][k] << "\t";
     461                output << endl;
     462        }
     463        output.close();
     464        return true;
    511465};
    512466
    513467// ======================================= CLASS EnergyMatrix =============================
     468
     469/** Create a trivial energy index mapping.
     470 * This just maps 1 to 1, 2 to 2 and so on for all fragments.
     471 * \return creation sucessful
     472 */
     473bool EnergyMatrix::ParseIndices()
     474{
     475        cout << "Parsing energy indices." << endl;
     476        Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "EnergyMatrix::ParseIndices: **Indices");
     477        for(int i=MatrixCounter+1;i--;) {
     478                Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "EnergyMatrix::ParseIndices: *Indices[]");
     479                for(int j=RowCounter[i];j--;)
     480                        Indices[i][j] = j;
     481        }
     482        return true;
     483};
    514484
    515485/** Sums the energy with each factor and put into last element of \a EnergyMatrix::Matrix.
    516486 * Sums over the "F"-terms in ANOVA decomposition.
    517  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     487 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
    518488 * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
    519489 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    … …  
    524494bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign)
    525495{
    526   // sum energy
    527   if (CorrectionFragments == NULL)
    528     for(int i=KeySet.FragmentsPerOrder[Order];i--;)
    529       for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
    530         for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
    531           Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
    532   else
    533     for(int i=KeySet.FragmentsPerOrder[Order];i--;)
    534       for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
    535         for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
    536           Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
    537   return true;
     496        // sum energy
     497        if (CorrectionFragments == NULL)
     498                for(int i=KeySet.FragmentsPerOrder[Order];i--;)
     499                        for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
     500                                for(int k=ColumnCounter;k--;)
     501                                        Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
     502        else
     503                for(int i=KeySet.FragmentsPerOrder[Order];i--;)
     504                        for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
     505                                for(int k=ColumnCounter;k--;)
     506                                        Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
     507        return true;
    538508};
    539509
    … …  
    545515 * \param skiplines number of inital columns to skip
    546516 * \return parsing successful
    547  */ 
     517 */
    548518bool EnergyMatrix::ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns)
    549519{
    550   char filename[1024];
    551   bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
    552 
    553   if (status) {
    554     // count maximum of columns
    555     RowCounter[MatrixCounter] = 0;
    556     ColumnCounter[MatrixCounter] = 0;
    557     for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
    558       if (RowCounter[j] > RowCounter[MatrixCounter])
    559         RowCounter[MatrixCounter] = RowCounter[j];
    560       if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
    561         ColumnCounter[MatrixCounter] = ColumnCounter[j];
    562     }
    563     // allocate last plus one matrix
    564     cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
    565     Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    566     for(int j=0;j<=RowCounter[MatrixCounter];j++)
    567       Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    568    
    569     // try independently to parse global energysuffix file if present
    570     strncpy(filename, name, 1023);
    571     strncat(filename, prefix, 1023-strlen(filename));
    572     strncat(filename, suffix.c_str(), 1023-strlen(filename));
    573     ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
    574   }
    575   return status;
     520        char filename[1024];
     521        bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
     522
     523        if (status) {
     524                // count maximum of columns
     525                RowCounter[MatrixCounter] = 0;
     526                for(int j=0; j < MatrixCounter;j++) // (energy matrix might be bigger than number of atoms in terms of rows)
     527                        if (RowCounter[j] > RowCounter[MatrixCounter])
     528                                RowCounter[MatrixCounter] = RowCounter[j];
     529                // allocate last plus one matrix
     530                cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
     531                Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     532                for(int j=0;j<=RowCounter[MatrixCounter];j++)
     533                        Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     534
     535                // try independently to parse global energysuffix file if present
     536                strncpy(filename, name, 1023);
     537                strncat(filename, prefix, 1023-strlen(filename));
     538                strncat(filename, suffix.c_str(), 1023-strlen(filename));
     539                ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
     540        }
     541        return status;
    576542};
    577543
    … …  
    580546/** Parsing force Indices of each fragment
    581547 * \param *name directory with \a ForcesFile
    582  * \return parsing successful 
    583  */
    584 bool ForceMatrix::ParseIndices(char *name) 
    585 {
    586   ifstream input;
    587   char *FragmentNumber = NULL;
    588   char filename[1023];
    589   stringstream line;
    590  
    591   cout << "Parsing force indices for " << MatrixCounter << " matrices." << endl;
    592   Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "ForceMatrix::ParseIndices: **Indices");
    593   line << name << FRAGMENTPREFIX << FORCESFILE;
    594   input.open(line.str().c_str(), ios::in);
    595   //cout << "Opening " << line.str() << " ... "  << input << endl;
    596   if (input == NULL) {
    597     cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
    598     return false;
    599   }
    600   for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
    601     // get the number of atoms for this fragment
    602     input.getline(filename, 1023);
    603     line.str(filename);
    604     // parse the values
    605     Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
    606     FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    607     //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
    608     Free((void **)&FragmentNumber, "ForceMatrix::ParseIndices: *FragmentNumber");
    609     for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
    610       line >> Indices[i][j];
    611       //cout << " " << Indices[i][j];
    612     }
    613     //cout << endl;
    614   }
    615   Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
    616   for(int j=RowCounter[MatrixCounter];j--;) {
    617     Indices[MatrixCounter][j] = j;
    618   }
    619   input.close();
    620   return true;
     548 * \return parsing successful
     549 */
     550bool ForceMatrix::ParseIndices(char *name)
     551{
     552        ifstream input;
     553        char *FragmentNumber = NULL;
     554        char filename[1023];
     555        stringstream line;
     556
     557        cout << "Parsing force indices for " << MatrixCounter << " matrices." << endl;
     558        Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "ForceMatrix::ParseIndices: **Indices");
     559        line << name << FRAGMENTPREFIX << FORCESFILE;
     560        input.open(line.str().c_str(), ios::in);
     561        //cout << "Opening " << line.str() << " ... "   << input << endl;
     562        if (input == NULL) {
     563                cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
     564                return false;
     565        }
     566        for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
     567                // get the number of atoms for this fragment
     568                input.getline(filename, 1023);
     569                line.str(filename);
     570                // parse the values
     571                Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
     572                FragmentNumber = FixedDigitNumber(MatrixCounter, i);
     573                //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
     574                Free((void **)&FragmentNumber, "ForceMatrix::ParseIndices: *FragmentNumber");
     575                for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
     576                        line >> Indices[i][j];
     577                        //cout << " " << Indices[i][j];
     578                }
     579                //cout << endl;
     580        }
     581        Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
     582        for(int j=RowCounter[MatrixCounter];j--;) {
     583                Indices[MatrixCounter][j] = j;
     584        }
     585        input.close();
     586        return true;
    621587};
    622588
    623589
    624590/** Sums the forces and puts into last element of \a ForceMatrix::Matrix.
    625  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     591 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
    626592 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    627593 * \param Order bond order
    628  *  \param sign +1 or -1
     594 *      \param sign +1 or -1
    629595 * \return true if summing was successful
    630596 */
    631597bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
    632598{
    633   int FragmentNr;
    634   // sum forces
    635   for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
    636     FragmentNr = KeySet.OrderSet[Order][i];
    637     for(int l=0;l<RowCounter[ FragmentNr ];l++) {
    638       int j = Indices[ FragmentNr ][l];
    639       if (j > RowCounter[MatrixCounter]) {
    640         cerr << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
    641         return false;
    642       }
    643       if (j != -1) {
    644         //if (j == 0) cout << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
    645         for(int k=2;k<ColumnCounter[MatrixCounter];k++)
    646           Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
    647       }
    648     }
    649   }
    650   return true;
     599        int FragmentNr;
     600        // sum forces
     601        for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
     602                FragmentNr = KeySet.OrderSet[Order][i];
     603                for(int l=0;l<RowCounter[ FragmentNr ];l++) {
     604                        int j = Indices[ FragmentNr ][l];
     605                        if (j > RowCounter[MatrixCounter]) {
     606                                cerr << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
     607                                return false;
     608                        }
     609                        if (j != -1) {
     610                                //if (j == 0) cout << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
     611                                for(int k=2;k<ColumnCounter;k++)
     612                                        Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
     613                        }
     614                }
     615        }
     616        return true;
    651617};
    652618
    … …  
    659625 * \param skiplines number of inital columns to skip
    660626 * \return parsing successful
    661  */ 
     627 */
    662628bool ForceMatrix::ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns)
    663629{
    664   char filename[1023];
    665   ifstream input;
    666   stringstream file;
    667   int nr;
    668   bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
    669 
    670   if (status) {
    671     // count number of atoms for last plus one matrix
    672     file << name << FRAGMENTPREFIX << KEYSETFILE;
    673     input.open(file.str().c_str(), ios::in);
    674     if (input == NULL) {
    675       cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
    676       return false;
    677     }
    678     RowCounter[MatrixCounter] = 0;
    679     while (!input.eof()) {
    680       input.getline(filename, 1023);
    681       stringstream zeile(filename);
    682       while (!zeile.eof()) {
    683         zeile >> nr;
    684         //cout << "Current index: " << nr << "." << endl;
    685         if (nr > RowCounter[MatrixCounter])
    686           RowCounter[MatrixCounter] = nr;
    687       }
    688     }
    689     RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
    690     input.close();
    691 
    692     ColumnCounter[MatrixCounter] = 0;
    693     for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
    694       if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
    695         ColumnCounter[MatrixCounter] = ColumnCounter[j];
    696     }
    697  
    698     // allocate last plus one matrix
    699     cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
    700     Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    701     for(int j=0;j<=RowCounter[MatrixCounter];j++)
    702       Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    703 
    704     // try independently to parse global forcesuffix file if present
    705     strncpy(filename, name, 1023);
    706     strncat(filename, prefix, 1023-strlen(filename));
    707     strncat(filename, suffix.c_str(), 1023-strlen(filename));
    708     ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
    709   }
    710  
    711 
    712   return status;
    713 };
    714 
    715 // ======================================= CLASS HessianMatrix =============================
    716 
    717 /** Parsing force Indices of each fragment
    718  * \param *name directory with \a ForcesFile
    719  * \return parsing successful
    720  */
    721 bool HessianMatrix::ParseIndices(char *name)
    722 {
    723   ifstream input;
    724   char *FragmentNumber = NULL;
    725   char filename[1023];
    726   stringstream line;
    727  
    728   cout << "Parsing hessian indices for " << MatrixCounter << " matrices." << endl;
    729   Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "HessianMatrix::ParseIndices: **Indices");
    730   line << name << FRAGMENTPREFIX << FORCESFILE;
    731   input.open(line.str().c_str(), ios::in);
    732   //cout << "Opening " << line.str() << " ... "  << input << endl;
    733   if (input == NULL) {
    734     cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
    735     return false;
    736   }
    737   for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
    738     // get the number of atoms for this fragment
    739     input.getline(filename, 1023);
    740     line.str(filename);
    741     // parse the values
    742     Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "HessianMatrix::ParseIndices: *Indices[]");
    743     FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    744     //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
    745     Free((void **)&FragmentNumber, "HessianMatrix::ParseIndices: *FragmentNumber");
    746     for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
    747       line >> Indices[i][j];
    748       //cout << " " << Indices[i][j];
    749     }
    750     //cout << endl;
    751   }
    752   Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "HessianMatrix::ParseIndices: *Indices[]");
    753   for(int j=RowCounter[MatrixCounter];j--;) {
    754     Indices[MatrixCounter][j] = j;
    755   }
    756   input.close();
    757   return true;
    758 };
    759 
    760 
    761 /** Sums the hessian entries and puts into last element of \a HessianMatrix::Matrix.
    762  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
    763  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    764  * \param Order bond order
    765  *  \param sign +1 or -1
    766  * \return true if summing was successful
    767  */
    768 bool HessianMatrix::SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
    769 {
    770   int FragmentNr;
    771   // sum forces
    772   for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
    773     FragmentNr = KeySet.OrderSet[Order][i];
    774     for(int l=0;l<RowCounter[ FragmentNr ];l++) {
    775       int j = Indices[ FragmentNr ][l];
    776       if (j > RowCounter[MatrixCounter]) {
    777         cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
    778         return false;
    779       }
    780       if (j != -1) {
    781         for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
    782           int k = Indices[ FragmentNr ][m];
    783           if (k > ColumnCounter[MatrixCounter]) {
    784             cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
    785             return false;
    786           }
    787           if (k != -1) {
    788             //cout << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
    789             Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
    790           }
    791         }
    792       }
    793     }
    794   }
    795   return true;
    796 };
    797 
    798 /** Constructor for class HessianMatrix.
    799  */
    800 HessianMatrix::HessianMatrix() : MatrixContainer()
    801 {
    802    IsSymmetric = true;
    803 }
    804 
    805 /** Sums the hessian entries with each factor and put into last element of \a ***Matrix.
    806  * Sums over "E"-terms to create the "F"-terms
    807  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
    808  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    809  * \param Order bond order
    810  * \return true if summing was successful
    811  */
    812 bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
    813 {
    814   // go through each order
    815   for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
    816     //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
    817     // then go per order through each suborder and pick together all the terms that contain this fragment
    818     for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
    819       for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
    820         if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
    821           //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
    822           // if the fragment's indices are all in the current fragment
    823           for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
    824             int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
    825             //cout << "Current row index is " << k << "/" << m << "." << endl;
    826             if (m != -1) { // if it's not an added hydrogen
    827               for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
    828                 //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
    829                 if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
    830                   m = l;
    831                   break; 
    832                 }
    833               }
    834               //cout << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
    835               if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
    836                 cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
    837                 return false;
    838               }
    839              
    840               for(int l=0;l<ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l++) {
    841                 int n = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][l];
    842                 //cout << "Current column index is " << l << "/" << n << "." << endl;
    843                 if (n != -1) { // if it's not an added hydrogen
    844                   for (int p=0;p<ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
    845                     //cout << "Comparing " << n << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
    846                     if (n == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p]) {
    847                       n = p;
    848                       break; 
    849                     }
    850                   }
    851                   //cout << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
    852                   if (n > ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
    853                     cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
    854                     return false;
    855                   }
    856                   if (Order == SubOrder) { // equal order is always copy from Energies
    857                     //cout << "Adding " << MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
    858                     Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
    859                   } else {
    860                     //cout << "Subtracting " << Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
    861                     Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
    862                   }
    863                 }
    864               }
    865             }
    866             //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
    867              //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
    868           }
    869         } else {
    870           //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
    871         }
    872       }
    873     }
    874    //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
    875   }
    876  
    877   return true;
    878 };
    879 
    880 /** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
    881  * \param *name directory with files
    882  * \param *prefix prefix of each matrix file
    883  * \param *suffix suffix of each matrix file
    884  * \param skiplines number of inital lines to skip
    885  * \param skiplines number of inital columns to skip
    886  * \return parsing successful
    887  */
    888 bool HessianMatrix::ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns)
    889 {
    890   char filename[1023];
    891   ifstream input;
    892   stringstream file;
    893   int nr;
    894   bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
    895 
    896   if (status) {
    897     // count number of atoms for last plus one matrix
    898     file << name << FRAGMENTPREFIX << KEYSETFILE;
    899     input.open(file.str().c_str(), ios::in);
    900     if (input == NULL) {
    901       cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
    902       return false;
    903     }
    904     RowCounter[MatrixCounter] = 0;
    905     ColumnCounter[MatrixCounter] = 0;
    906     while (!input.eof()) {
    907       input.getline(filename, 1023);
    908       stringstream zeile(filename);
    909       while (!zeile.eof()) {
    910         zeile >> nr;
    911         //cout << "Current index: " << nr << "." << endl;
    912         if (nr > RowCounter[MatrixCounter]) {
    913           RowCounter[MatrixCounter] = nr;
    914           ColumnCounter[MatrixCounter] = nr;
    915         }
    916       }
    917     }
    918     RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
    919     ColumnCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
    920     input.close();
    921  
    922     // allocate last plus one matrix
    923     cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
    924     Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    925     for(int j=0;j<=RowCounter[MatrixCounter];j++)
    926       Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    927 
    928     // try independently to parse global forcesuffix file if present
    929     strncpy(filename, name, 1023);
    930     strncat(filename, prefix, 1023-strlen(filename));
    931     strncat(filename, suffix.c_str(), 1023-strlen(filename));
    932     ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
    933   }
    934  
    935 
    936   return status;
     630        char filename[1023];
     631        ifstream input;
     632        stringstream file;
     633        int nr;
     634        bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
     635
     636        if (status) {
     637                // count number of atoms for last plus one matrix
     638                file << name << FRAGMENTPREFIX << KEYSETFILE;
     639                input.open(file.str().c_str(), ios::in);
     640                if (input == NULL) {
     641                        cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
     642                        return false;
     643                }
     644                RowCounter[MatrixCounter] = 0;
     645                while (!input.eof()) {
     646                        input.getline(filename, 1023);
     647                        stringstream zeile(filename);
     648                        while (!zeile.eof()) {
     649                                zeile >> nr;
     650                                //cout << "Current index: " << nr << "." << endl;
     651                                if (nr > RowCounter[MatrixCounter])
     652                                        RowCounter[MatrixCounter] = nr;
     653                        }
     654                }
     655                RowCounter[MatrixCounter]++;            // nr start at 0, count starts at 1
     656                input.close();
     657
     658                // allocate last plus one matrix
     659                cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
     660                Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     661                for(int j=0;j<=RowCounter[MatrixCounter];j++)
     662                        Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     663
     664                // try independently to parse global forcesuffix file if present
     665                strncpy(filename, name, 1023);
     666                strncat(filename, prefix, 1023-strlen(filename));
     667                strncat(filename, suffix.c_str(), 1023-strlen(filename));
     668                ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
     669        }
     670
     671
     672        return status;
    937673};
    938674
    … …  
    942678 */
    943679KeySetsContainer::KeySetsContainer() {
    944   KeySets = NULL;
    945   AtomCounter = NULL;
    946   FragmentCounter = 0;
    947   Order = 0;
    948   FragmentsPerOrder = 0;
    949   OrderSet = NULL;
     680        KeySets = NULL;
     681        AtomCounter = NULL;
     682        FragmentCounter = 0;
     683        Order = 0;
     684        FragmentsPerOrder = 0;
     685        OrderSet = NULL;
    950686};
    951687
    … …  
    953689 */
    954690KeySetsContainer::~KeySetsContainer() {
    955   for(int i=FragmentCounter;i--;)
    956     Free((void **)&KeySets[i], "KeySetsContainer::~KeySetsContainer: *KeySets[]");
    957   for(int i=Order;i--;)
    958     Free((void **)&OrderSet[i], "KeySetsContainer::~KeySetsContainer: *OrderSet[]");
    959   Free((void **)&KeySets, "KeySetsContainer::~KeySetsContainer: **KeySets");
    960   Free((void **)&OrderSet, "KeySetsContainer::~KeySetsContainer: **OrderSet");
    961   Free((void **)&AtomCounter, "KeySetsContainer::~KeySetsContainer: *AtomCounter");
    962   Free((void **)&FragmentsPerOrder, "KeySetsContainer::~KeySetsContainer: *FragmentsPerOrder");
     691        for(int i=FragmentCounter;i--;)
     692                Free((void **)&KeySets[i], "KeySetsContainer::~KeySetsContainer: *KeySets[]");
     693        for(int i=Order;i--;)
     694                Free((void **)&OrderSet[i], "KeySetsContainer::~KeySetsContainer: *OrderSet[]");
     695        Free((void **)&KeySets, "KeySetsContainer::~KeySetsContainer: **KeySets");
     696        Free((void **)&OrderSet, "KeySetsContainer::~KeySetsContainer: **OrderSet");
     697        Free((void **)&AtomCounter, "KeySetsContainer::~KeySetsContainer: *AtomCounter");
     698        Free((void **)&FragmentsPerOrder, "KeySetsContainer::~KeySetsContainer: *FragmentsPerOrder");
    963699};
    964700
    … …  
    970706 */
    971707bool KeySetsContainer::ParseKeySets(const char *name, const int *ACounter, const int FCounter) {
    972   ifstream input;
    973   char *FragmentNumber = NULL;
    974   stringstream file;
    975   char filename[1023];
    976  
    977   FragmentCounter = FCounter;
    978   cout << "Parsing key sets." << endl;
    979   KeySets = (int **) Malloc(sizeof(int *)*FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
    980   for(int i=FragmentCounter;i--;)
    981     KeySets[i] = NULL;
    982   file << name << FRAGMENTPREFIX << KEYSETFILE;
    983   input.open(file.str().c_str(), ios::in);
    984   if (input == NULL) {
    985     cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
    986     return false;
    987   }
    988  
    989   AtomCounter = (int *) Malloc(sizeof(int)*FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
    990   for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
    991     stringstream line;
    992     AtomCounter[i] = ACounter[i];
    993     // parse the values
    994     KeySets[i] = (int *) Malloc(sizeof(int)*AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
    995     for(int j=AtomCounter[i];j--;)
    996       KeySets[i][j] = -1;
    997     FragmentNumber = FixedDigitNumber(FragmentCounter, i);
    998     //cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
    999     Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
    1000     input.getline(filename, 1023);
    1001     line.str(filename);
    1002     for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
    1003       line >> KeySets[i][j];
    1004       //cout << " " << KeySets[i][j];
    1005     }
    1006     //cout << endl;
    1007   }
    1008   input.close();
    1009   return true;
     708        ifstream input;
     709        char *FragmentNumber = NULL;
     710        stringstream file;
     711        char filename[1023];
     712
     713        FragmentCounter = FCounter;
     714        cout << "Parsing key sets." << endl;
     715        KeySets = (int **) Malloc(sizeof(int *)*FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
     716        for(int i=FragmentCounter;i--;)
     717                KeySets[i] = NULL;
     718        file << name << FRAGMENTPREFIX << KEYSETFILE;
     719        input.open(file.str().c_str(), ios::in);
     720        if (input == NULL) {
     721                cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
     722                return false;
     723        }
     724
     725        AtomCounter = (int *) Malloc(sizeof(int)*FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
     726        for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
     727                stringstream line;
     728                AtomCounter[i] = ACounter[i];
     729                // parse the values
     730                KeySets[i] = (int *) Malloc(sizeof(int)*AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
     731                for(int j=AtomCounter[i];j--;)
     732                        KeySets[i][j] = -1;
     733                FragmentNumber = FixedDigitNumber(FragmentCounter, i);
     734                //cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
     735                Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
     736                input.getline(filename, 1023);
     737                line.str(filename);
     738                for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
     739                        line >> KeySets[i][j];
     740                        //cout << " " << KeySets[i][j];
     741                }
     742                //cout << endl;
     743        }
     744        input.close();
     745        return true;
    1010746};
    1011747
    … …  
    1015751bool KeySetsContainer::ParseManyBodyTerms()
    1016752{
    1017   int Counter;
    1018  
    1019   cout << "Creating Fragment terms." << endl;
    1020   // scan through all to determine maximum order
    1021   Order=0;
    1022   for(int i=FragmentCounter;i--;) {
    1023     Counter=0;
    1024     for(int j=AtomCounter[i];j--;)
    1025       if (KeySets[i][j] != -1)
    1026         Counter++;
    1027     if (Counter > Order)
    1028       Order = Counter;
    1029   }
    1030   cout << "Found Order is " << Order << "." << endl;
    1031  
    1032   // scan through all to determine fragments per order
    1033   FragmentsPerOrder = (int *) Malloc(sizeof(int)*Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
    1034   for(int i=Order;i--;)
    1035     FragmentsPerOrder[i] = 0;
    1036   for(int i=FragmentCounter;i--;) {
    1037     Counter=0;
    1038     for(int j=AtomCounter[i];j--;)
    1039       if (KeySets[i][j] != -1)
    1040         Counter++;
    1041     FragmentsPerOrder[Counter-1]++;
    1042   }
    1043   for(int i=0;i<Order;i++)
    1044     cout << "Found No. of Fragments of Order " << i+1 << " is " << FragmentsPerOrder[i] << "." << endl;
    1045    
    1046   // scan through all to gather indices to each order set
    1047   OrderSet = (int **) Malloc(sizeof(int *)*Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
    1048   for(int i=Order;i--;)
    1049     OrderSet[i] = (int *) Malloc(sizeof(int)*FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
    1050   for(int i=Order;i--;)
    1051     FragmentsPerOrder[i] = 0;
    1052   for(int i=FragmentCounter;i--;) {
    1053     Counter=0;
    1054     for(int j=AtomCounter[i];j--;)
    1055       if (KeySets[i][j] != -1)
    1056         Counter++;
    1057     OrderSet[Counter-1][FragmentsPerOrder[Counter-1]] = i;
    1058     FragmentsPerOrder[Counter-1]++;
    1059   }
    1060   cout << "Printing OrderSet." << endl;
    1061   for(int i=0;i<Order;i++) {
    1062     for (int j=0;j<FragmentsPerOrder[i];j++) {
    1063       cout << " " << OrderSet[i][j];
    1064     }
    1065     cout << endl;
    1066   }
    1067   cout << endl;
    1068  
    1069    
    1070   return true;
     753        int Counter;
     754
     755        cout << "Creating Fragment terms." << endl;
     756        // scan through all to determine maximum order
     757        Order=0;
     758        for(int i=FragmentCounter;i--;) {
     759                Counter=0;
     760                for(int j=AtomCounter[i];j--;)
     761                        if (KeySets[i][j] != -1)
     762                                Counter++;
     763                if (Counter > Order)
     764                        Order = Counter;
     765        }
     766        cout << "Found Order is " << Order << "." << endl;
     767
     768        // scan through all to determine fragments per order
     769        FragmentsPerOrder = (int *) Malloc(sizeof(int)*Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
     770        for(int i=Order;i--;)
     771                FragmentsPerOrder[i] = 0;
     772        for(int i=FragmentCounter;i--;) {
     773                Counter=0;
     774                for(int j=AtomCounter[i];j--;)
     775                        if (KeySets[i][j] != -1)
     776                                Counter++;
     777                FragmentsPerOrder[Counter-1]++;
     778        }
     779        for(int i=0;i<Order;i++)
     780                cout << "Found No. of Fragments of Order " << i+1 << " is " << FragmentsPerOrder[i] << "." << endl;
     781
     782        // scan through all to gather indices to each order set
     783        OrderSet = (int **) Malloc(sizeof(int *)*Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
     784        for(int i=Order;i--;)
     785                OrderSet[i] = (int *) Malloc(sizeof(int)*FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
     786        for(int i=Order;i--;)
     787                FragmentsPerOrder[i] = 0;
     788        for(int i=FragmentCounter;i--;) {
     789                Counter=0;
     790                for(int j=AtomCounter[i];j--;)
     791                        if (KeySets[i][j] != -1)
     792                                Counter++;
     793                OrderSet[Counter-1][FragmentsPerOrder[Counter-1]] = i;
     794                FragmentsPerOrder[Counter-1]++;
     795        }
     796        cout << "Printing OrderSet." << endl;
     797        for(int i=0;i<Order;i++) {
     798                for (int j=0;j<FragmentsPerOrder[i];j++) {
     799                        cout << " " << OrderSet[i][j];
     800                }
     801                cout << endl;
     802        }
     803        cout << endl;
     804
     805
     806        return true;
    1071807};
    1072808
    … …  
    1078814bool KeySetsContainer::Contains(const int GreaterSet, const int SmallerSet)
    1079815{
    1080   bool result = true;
    1081   bool intermediate;
    1082   if ((GreaterSet < 0) || (SmallerSet < 0) || (GreaterSet > FragmentCounter) || (SmallerSet > FragmentCounter)) // index out of bounds
    1083     return false;
    1084   for(int i=AtomCounter[SmallerSet];i--;) {
    1085     intermediate = false;
    1086     for (int j=AtomCounter[GreaterSet];j--;)
    1087       intermediate = (intermediate || ((KeySets[SmallerSet][i] == KeySets[GreaterSet][j]) || (KeySets[SmallerSet][i] == -1)));
    1088     result = result && intermediate;
    1089   }
    1090 
    1091   return result;
     816        bool result = true;
     817        bool intermediate;
     818        if ((GreaterSet < 0) || (SmallerSet < 0) || (GreaterSet > FragmentCounter) || (SmallerSet > FragmentCounter)) // index out of bounds
     819                return false;
     820        for(int i=AtomCounter[SmallerSet];i--;) {
     821                intermediate = false;
     822                for (int j=AtomCounter[GreaterSet];j--;)
     823                        intermediate = (intermediate || ((KeySets[SmallerSet][i] == KeySets[GreaterSet][j]) || (KeySets[SmallerSet][i] == -1)));
     824                result = result && intermediate;
     825        }
     826
     827        return result;
    1092828};
    1093829
Note: See TracChangeset for help on using the changeset viewer.