Ignore:
Timestamp:
Feb 9, 2009, 5:24:10 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
451d7a
Parents:
4aef8a
git-author:
Frederik Heber <heber@…> (02/09/09 15:55:37)
git-committer:
Frederik Heber <heber@…> (02/09/09 17:24:10)
Message:

Merge branch 'ConcaveHull' of ../espack2 into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/linkedcell.hpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp
util/src/NanoCreator.c

Basically, this resulted from a lot of conversions two from spaces to one tab, which is my standard indentation. The mess was caused by eclipse auto-indenting. And in espack2:ConcaveHull was the new stuff, so all from ConcaveHull was replaced in case of doubt.
Additionally, vector had ofstream << operator instead ostream << ...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/moleculelist.cpp

    • Property mode changed from 100644 to 100755
    r4aef8a re08f45  
    2121MoleculeListClass::MoleculeListClass(int NumMolecules, int NumAtoms)
    2222{
    23   ListOfMolecules = (molecule **) Malloc(sizeof(molecule *)*NumMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
    24   for (int i=NumMolecules;i--;)
    25     ListOfMolecules[i] = NULL;
    26   NumberOfMolecules = NumMolecules;
    27   NumberOfTopAtoms = NumAtoms;
     23        ListOfMolecules = (molecule **) Malloc(sizeof(molecule *)*NumMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
     24        for (int i=NumMolecules;i--;)
     25                ListOfMolecules[i] = NULL;
     26        NumberOfMolecules = NumMolecules;
     27        NumberOfTopAtoms = NumAtoms;
    2828};
    2929
    … …  
    3333MoleculeListClass::~MoleculeListClass()
    3434{
    35   cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
    36   for (int i=NumberOfMolecules;i--;) {
    37     if (ListOfMolecules[i] != NULL) { // if NULL don't free
    38       cout << Verbose(4) << "ListOfMolecules: Freeing " << ListOfMolecules[i] << "." << endl;
    39       delete(ListOfMolecules[i]);
    40       ListOfMolecules[i] = NULL;
    41     }
    42   }
    43   cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
    44   Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
     35        cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
     36        for (int i=NumberOfMolecules;i--;) {
     37                if (ListOfMolecules[i] != NULL) { // if NULL don't free
     38                        cout << Verbose(4) << "ListOfMolecules: Freeing " << ListOfMolecules[i] << "." << endl;
     39                        delete(ListOfMolecules[i]);
     40                        ListOfMolecules[i] = NULL;
     41                }
     42        }
     43        cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
     44        Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
    4545};
    4646
    … …  
    5252int MolCompare(const void *a, const void *b)
    5353{
    54   int *aList = NULL, *bList = NULL;
    55   int Count, Counter, aCounter, bCounter;
    56   int flag;
    57   atom *aWalker = NULL;
    58   atom *bWalker = NULL;
    59  
    60   // sort each atom list and put the numbers into a list, then go through
    61   //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
    62   if ( (**(molecule **)a).AtomCount < (**(molecule **)b).AtomCount ) {
    63     return -1;
    64   } else { if ((**(molecule **)a).AtomCount > (**(molecule **)b).AtomCount)
    65     return +1;
    66     else {
    67       Count = (**(molecule **)a).AtomCount;
    68       aList = new int[Count];
    69       bList = new int[Count];
    70  
    71       // fill the lists
    72       aWalker = (**(molecule **)a).start;
    73       bWalker = (**(molecule **)b).start;
    74       Counter = 0;
    75       aCounter = 0;
    76       bCounter = 0;
    77       while ((aWalker->next != (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
    78         aWalker = aWalker->next;
    79         bWalker = bWalker->next;
    80         if (aWalker->GetTrueFather() == NULL)
    81           aList[Counter] = Count + (aCounter++);
    82         else
    83           aList[Counter] = aWalker->GetTrueFather()->nr;
    84         if (bWalker->GetTrueFather() == NULL)
    85           bList[Counter] = Count + (bCounter++);
    86         else
    87           bList[Counter] = bWalker->GetTrueFather()->nr;
    88         Counter++;
    89       }
    90       // check if AtomCount was for real
    91       flag = 0;
    92       if ((aWalker->next == (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
    93         flag = -1;
    94       } else {
    95         if ((aWalker->next != (**(molecule **)a).end) && (bWalker->next == (**(molecule **)b).end))
    96           flag = 1;
    97       }
    98       if (flag == 0) {
    99         // sort the lists
    100         gsl_heapsort(aList,Count, sizeof(int), CompareDoubles);
    101         gsl_heapsort(bList,Count, sizeof(int), CompareDoubles);
    102         // compare the lists
    103        
    104         flag = 0;
    105         for(int i=0;i<Count;i++) {
    106           if (aList[i] < bList[i]) {
    107             flag = -1;
    108           } else {
    109             if (aList[i] > bList[i])
    110               flag = 1;
    111           }
    112           if (flag != 0)
    113             break;
    114         }
    115       }
    116       delete[](aList);
    117       delete[](bList);
    118       return flag;
    119     }
    120   }
    121   return  -1;
     54        int *aList = NULL, *bList = NULL;
     55        int Count, Counter, aCounter, bCounter;
     56        int flag;
     57        atom *aWalker = NULL;
     58        atom *bWalker = NULL;
     59       
     60        // sort each atom list and put the numbers into a list, then go through
     61        //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
     62        if ( (**(molecule **)a).AtomCount < (**(molecule **)b).AtomCount ) {
     63                return -1;
     64        } else { if ((**(molecule **)a).AtomCount > (**(molecule **)b).AtomCount)
     65                return +1;
     66                else {
     67                        Count = (**(molecule **)a).AtomCount;
     68                        aList = new int[Count];
     69                        bList = new int[Count];
     70       
     71                        // fill the lists
     72                        aWalker = (**(molecule **)a).start;
     73                        bWalker = (**(molecule **)b).start;
     74                        Counter = 0;
     75                        aCounter = 0;
     76                        bCounter = 0;
     77                        while ((aWalker->next != (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
     78                                aWalker = aWalker->next;
     79                                bWalker = bWalker->next;
     80                                if (aWalker->GetTrueFather() == NULL)
     81                                        aList[Counter] = Count + (aCounter++);
     82                                else
     83                                        aList[Counter] = aWalker->GetTrueFather()->nr;
     84                                if (bWalker->GetTrueFather() == NULL)
     85                                        bList[Counter] = Count + (bCounter++);
     86                                else
     87                                        bList[Counter] = bWalker->GetTrueFather()->nr;
     88                                Counter++;
     89                        }
     90                        // check if AtomCount was for real
     91                        flag = 0;
     92                        if ((aWalker->next == (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
     93                                flag = -1;
     94                        } else {
     95                                if ((aWalker->next != (**(molecule **)a).end) && (bWalker->next == (**(molecule **)b).end))
     96                                        flag = 1;
     97                        }
     98                        if (flag == 0) {
     99                                // sort the lists
     100                                gsl_heapsort(aList,Count, sizeof(int), CompareDoubles);
     101                                gsl_heapsort(bList,Count, sizeof(int), CompareDoubles);
     102                                // compare the lists
     103                               
     104                                flag = 0;
     105                                for(int i=0;i<Count;i++) {
     106                                        if (aList[i] < bList[i]) {
     107                                                flag = -1;
     108                                        } else {
     109                                                if (aList[i] > bList[i])
     110                                                        flag = 1;
     111                                        }
     112                                        if (flag != 0)
     113                                                break;
     114                                }
     115                        }
     116                        delete[](aList);
     117                        delete[](bList);
     118                        return flag;
     119                }
     120        }
     121        return  -1;
    122122};
    123123
    … …  
    127127void MoleculeListClass::Output(ofstream *out)
    128128{
    129   *out<< Verbose(1) << "MoleculeList: ";
    130   for (int i=0;i<NumberOfMolecules;i++)
    131     *out << ListOfMolecules[i] << "\t";
    132   *out << endl;
     129        *out<< Verbose(1) << "MoleculeList: ";
     130        for (int i=0;i<NumberOfMolecules;i++)
     131                *out << ListOfMolecules[i] << "\t";
     132        *out << endl;
    133133};
    134134
    … …  
    142142bool MoleculeListClass::AddHydrogenCorrection(ofstream *out, char *path)
    143143{
    144   atom *Walker = NULL;
    145   atom *Runner = NULL;
    146   double ***FitConstant = NULL, **correction = NULL;
    147   int a,b;
    148   ofstream output;
    149   ifstream input;
    150   string line;
    151   stringstream zeile;
    152   double distance;
    153   char ParsedLine[1023];
    154   double tmp;
    155   char *FragmentNumber = NULL;
    156 
    157   cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
    158   // 0. parse in fit constant files that should have the same dimension as the final energy files
    159   // 0a. find dimension of matrices with constants
    160   line = path;
    161   line.append("/");
    162   line += FRAGMENTPREFIX;
    163   line += "1";
    164   line += FITCONSTANTSUFFIX;
    165   input.open(line.c_str());
    166   if (input == NULL) {
    167     cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
    168     return false;
    169   }
    170   a=0;
    171   b=-1; // we overcount by one
    172   while (!input.eof()) {
    173     input.getline(ParsedLine, 1023);
    174     zeile.str(ParsedLine);
    175     int i=0;
    176     while (!zeile.eof()) {
    177       zeile >> distance;
    178       i++;
    179     }
    180     if (i > a)
    181       a = i;
    182     b++;
    183   }
    184   cout << "I recognized " << a << " columns and " << b << " rows, ";
    185   input.close();
    186  
    187   // 0b. allocate memory for constants
    188   FitConstant = (double ***) Malloc(sizeof(double **)*3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
    189   for (int k=0;k<3;k++) {
    190     FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
    191     for (int i=a;i--;) {
    192       FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
    193     }
    194   }
    195   // 0c. parse in constants
    196   for (int i=0;i<3;i++) {
    197     line = path;
    198     line.append("/");
    199     line += FRAGMENTPREFIX;
    200     sprintf(ParsedLine, "%d", i+1);
    201     line += ParsedLine;
    202     line += FITCONSTANTSUFFIX;
    203     input.open(line.c_str());
    204     if (input == NULL) {
    205       cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
    206       return false;
    207     }
    208     int k = 0,l;
    209     while ((!input.eof()) && (k < b)) {
    210       input.getline(ParsedLine, 1023);
    211       //cout << "Current Line: " << ParsedLine << endl;
    212       zeile.str(ParsedLine);
    213       zeile.clear();
    214       l = 0;
    215       while ((!zeile.eof()) && (l < a)) {
    216         zeile >> FitConstant[i][l][k];
    217         //cout << FitConstant[i][l][k] << "\t";
    218         l++;
    219       }
    220       //cout << endl;
    221       k++;
    222     }
    223     input.close();
    224   }
    225   for(int k=0;k<3;k++) {
    226     cout << "Constants " << k << ":" << endl;
    227     for (int j=0;j<b;j++) {
    228       for (int i=0;i<a;i++) {
    229         cout << FitConstant[k][i][j] << "\t";
    230       }
    231       cout << endl;
    232     }
    233     cout << endl;
    234   }
    235  
    236   // 0d. allocate final correction matrix
    237   correction = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **correction");
    238   for (int i=a;i--;)
    239     correction[i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
    240        
    241   // 1a. go through every molecule in the list
    242   for(int i=NumberOfMolecules;i--;) {
    243     // 1b. zero final correction matrix
    244     for (int k=a;k--;)
    245       for (int j=b;j--;)
    246         correction[k][j] = 0.;
    247     // 2. take every hydrogen that is a saturated one
    248     Walker = ListOfMolecules[i]->start;
    249     while (Walker->next != ListOfMolecules[i]->end) {
    250       Walker = Walker->next;
    251       //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
    252       if ((Walker->type->Z == 1) && ((Walker->father == NULL) || (Walker->father->type->Z != 1))) { // if it's a hydrogen
    253         Runner = ListOfMolecules[i]->start;
    254         while (Runner->next != ListOfMolecules[i]->end) {
    255           Runner = Runner->next;
    256           //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
    257           // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
    258           if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) {  // (hydrogens have only one bonding partner!)
    259             // 4. evaluate the morse potential for each matrix component and add up
    260             distance = sqrt(Runner->x.Distance(&Walker->x));
    261             //cout << "Fragment " << i << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
    262             for(int k=0;k<a;k++) {
    263               for (int j=0;j<b;j++) {
    264                 switch(k) {
    265                 case 1:
    266                 case 7:
    267                 case 11:
    268                   tmp = pow(FitConstant[0][k][j] * ( 1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]) ) ),2);
    269                   break;
    270                 default:
    271                   tmp = FitConstant[0][k][j] * pow( distance,  FitConstant[1][k][j]) + FitConstant[2][k][j];
    272                 };
    273                 correction[k][j] -= tmp;    // ground state is actually lower (disturbed by additional interaction)
    274                 //cout << tmp << "\t";
    275               }
    276               //cout << endl;
    277             }
    278             //cout << endl;
    279           }
    280         }
    281       }
    282     }
    283     // 5. write final matrix to file
    284     line = path;
    285     line.append("/");
    286     line += FRAGMENTPREFIX;
    287     FragmentNumber = FixedDigitNumber(NumberOfMolecules, i);
    288     line += FragmentNumber;
    289     delete(FragmentNumber);
    290     line += HCORRECTIONSUFFIX;
    291     output.open(line.c_str());
    292     output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
    293     for (int j=0;j<b;j++) {
    294       for(int i=0;i<a;i++)
    295         output << correction[i][j] << "\t";
    296       output << endl;
    297     }
    298     output.close();
    299   }
    300   line = path;
    301   line.append("/");
    302   line += HCORRECTIONSUFFIX;
    303   output.open(line.c_str());
    304   output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
    305   for (int j=0;j<b;j++) {
    306     for(int i=0;i<a;i++)
    307       output << 0 << "\t";
    308     output << endl;
    309   }
    310   output.close();
    311   // 6. free memory of parsed matrices
    312   FitConstant = (double ***) Malloc(sizeof(double **)*a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
    313   for (int k=0;k<3;k++) {
    314     FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
    315     for (int i=a;i--;) {
    316       FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
    317     }
    318   }
    319   cout << "done." << endl;
    320   return true;
     144        atom *Walker = NULL;
     145        atom *Runner = NULL;
     146        double ***FitConstant = NULL, **correction = NULL;
     147        int a,b;
     148        ofstream output;
     149        ifstream input;
     150        string line;
     151        stringstream zeile;
     152        double distance;
     153        char ParsedLine[1023];
     154        double tmp;
     155        char *FragmentNumber = NULL;
     156
     157        cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
     158        // 0. parse in fit constant files that should have the same dimension as the final energy files
     159        // 0a. find dimension of matrices with constants
     160        line = path;
     161        line.append("/");
     162        line += FRAGMENTPREFIX;
     163        line += "1";
     164        line += FITCONSTANTSUFFIX;
     165        input.open(line.c_str());
     166        if (input == NULL) {
     167                cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
     168                return false;
     169        }
     170        a=0;
     171        b=-1; // we overcount by one
     172        while (!input.eof()) {
     173                input.getline(ParsedLine, 1023);
     174                zeile.str(ParsedLine);
     175                int i=0;
     176                while (!zeile.eof()) {
     177                        zeile >> distance;
     178                        i++;
     179                }
     180                if (i > a)
     181                        a = i;
     182                b++;
     183        }
     184        cout << "I recognized " << a << " columns and " << b << " rows, ";
     185        input.close();
     186       
     187        // 0b. allocate memory for constants
     188        FitConstant = (double ***) Malloc(sizeof(double **)*3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
     189        for (int k=0;k<3;k++) {
     190                FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
     191                for (int i=a;i--;) {
     192                        FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
     193                }
     194        }
     195        // 0c. parse in constants
     196        for (int i=0;i<3;i++) {
     197                line = path;
     198                line.append("/");
     199                line += FRAGMENTPREFIX;
     200                sprintf(ParsedLine, "%d", i+1);
     201                line += ParsedLine;
     202                line += FITCONSTANTSUFFIX;
     203                input.open(line.c_str());
     204                if (input == NULL) {
     205                        cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
     206                        return false;
     207                }
     208                int k = 0,l;
     209                while ((!input.eof()) && (k < b)) {
     210                        input.getline(ParsedLine, 1023);
     211                        //cout << "Current Line: " << ParsedLine << endl;
     212                        zeile.str(ParsedLine);
     213                        zeile.clear();
     214                        l = 0;
     215                        while ((!zeile.eof()) && (l < a)) {
     216                                zeile >> FitConstant[i][l][k];
     217                                //cout << FitConstant[i][l][k] << "\t";
     218                                l++;
     219                        }
     220                        //cout << endl;
     221                        k++;
     222                }
     223                input.close();
     224        }
     225        for(int k=0;k<3;k++) {
     226                cout << "Constants " << k << ":" << endl;
     227                for (int j=0;j<b;j++) {
     228                        for (int i=0;i<a;i++) {
     229                                cout << FitConstant[k][i][j] << "\t";
     230                        }
     231                        cout << endl;
     232                }
     233                cout << endl;
     234        }
     235       
     236        // 0d. allocate final correction matrix
     237        correction = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **correction");
     238        for (int i=a;i--;)
     239                correction[i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
     240                               
     241        // 1a. go through every molecule in the list
     242        for(int i=NumberOfMolecules;i--;) {
     243                // 1b. zero final correction matrix
     244                for (int k=a;k--;)
     245                        for (int j=b;j--;)
     246                                correction[k][j] = 0.;
     247                // 2. take every hydrogen that is a saturated one
     248                Walker = ListOfMolecules[i]->start;
     249                while (Walker->next != ListOfMolecules[i]->end) {
     250                        Walker = Walker->next;
     251                        //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
     252                        if ((Walker->type->Z == 1) && ((Walker->father == NULL) || (Walker->father->type->Z != 1))) { // if it's a hydrogen
     253                                Runner = ListOfMolecules[i]->start;
     254                                while (Runner->next != ListOfMolecules[i]->end) {
     255                                        Runner = Runner->next;
     256                                        //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
     257                                        // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
     258                                        if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) {      // (hydrogens have only one bonding partner!)
     259                                                // 4. evaluate the morse potential for each matrix component and add up
     260                                                distance = Runner->x.Distance(&Walker->x);
     261                                                //cout << "Fragment " << i << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
     262                                                for(int k=0;k<a;k++) {
     263                                                        for (int j=0;j<b;j++) {
     264                                                                switch(k) {
     265                                                                case 1:
     266                                                                case 7:
     267                                                                case 11:
     268                                                                        tmp = pow(FitConstant[0][k][j] * ( 1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]) ) ),2);
     269                                                                        break;
     270                                                                default:
     271                                                                        tmp = FitConstant[0][k][j] * pow( distance,     FitConstant[1][k][j]) + FitConstant[2][k][j];
     272                                                                };
     273                                                                correction[k][j] -= tmp;                // ground state is actually lower (disturbed by additional interaction)
     274                                                                //cout << tmp << "\t";
     275                                                        }
     276                                                        //cout << endl;
     277                                                }
     278                                                //cout << endl;
     279                                        }
     280                                }
     281                        }
     282                }
     283                // 5. write final matrix to file
     284                line = path;
     285                line.append("/");
     286                line += FRAGMENTPREFIX;
     287                FragmentNumber = FixedDigitNumber(NumberOfMolecules, i);
     288                line += FragmentNumber;
     289                delete(FragmentNumber);
     290                line += HCORRECTIONSUFFIX;
     291                output.open(line.c_str());
     292                output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
     293                for (int j=0;j<b;j++) {
     294                        for(int i=0;i<a;i++)
     295                                output << correction[i][j] << "\t";
     296                        output << endl;
     297                }
     298                output.close();
     299        }
     300        line = path;
     301        line.append("/");
     302        line += HCORRECTIONSUFFIX;
     303        output.open(line.c_str());
     304        output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
     305        for (int j=0;j<b;j++) {
     306                for(int i=0;i<a;i++)
     307                        output << 0 << "\t";
     308                output << endl;
     309        }
     310        output.close();
     311        // 6. free memory of parsed matrices
     312        FitConstant = (double ***) Malloc(sizeof(double **)*a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
     313        for (int k=0;k<3;k++) {
     314                FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
     315                for (int i=a;i--;) {
     316                        FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
     317                }
     318        }
     319        cout << "done." << endl;
     320        return true;
    321321};
    322322
    … …  
    329329bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path, int *SortIndex)
    330330{
    331   bool status = true;
    332   ofstream ForcesFile;
    333   stringstream line;
    334   atom *Walker = NULL;
    335   element *runner = NULL;
    336 
    337   // open file for the force factors
    338   *out << Verbose(1) << "Saving  force factors ... ";
    339   line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
    340   ForcesFile.open(line.str().c_str(), ios::out);
    341   if (ForcesFile != NULL) {
    342     //cout << Verbose(1) << "Final AtomicForcesList: ";
    343     //output << prefix << "Forces" << endl;
    344     for(int j=0;j<NumberOfMolecules;j++) {
    345       //if (TEList[j] != 0) {
    346       runner = ListOfMolecules[j]->elemente->start;
    347       while (runner->next != ListOfMolecules[j]->elemente->end) { // go through every element
    348         runner = runner->next;
    349         if (ListOfMolecules[j]->ElementsInMolecule[runner->Z]) { // if this element got atoms
    350           Walker = ListOfMolecules[j]->start;
    351           while (Walker->next != ListOfMolecules[j]->end) { // go through every atom of this element
    352             Walker = Walker->next;
    353             if (Walker->type->Z == runner->Z) {
    354               if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
    355                 //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
    356                 ForcesFile <<  SortIndex[Walker->GetTrueFather()->nr] << "\t";
    357                 } else  // otherwise a -1 to indicate an added saturation hydrogen
    358                   ForcesFile << "-1\t";
    359               }
    360             }
    361           }
    362       }
    363       ForcesFile << endl;
    364     }
    365     ForcesFile.close();
    366     *out << Verbose(1) << "done." << endl;
    367   } else {
    368     status = false;
    369     *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
    370   }
    371   ForcesFile.close();
    372  
    373   return status;
     331        bool status = true;
     332        ofstream ForcesFile;
     333        stringstream line;
     334        atom *Walker = NULL;
     335        element *runner = NULL;
     336
     337        // open file for the force factors
     338        *out << Verbose(1) << "Saving   force factors ... ";
     339        line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
     340        ForcesFile.open(line.str().c_str(), ios::out);
     341        if (ForcesFile != NULL) {
     342                //cout << Verbose(1) << "Final AtomicForcesList: ";
     343                //output << prefix << "Forces" << endl;
     344                for(int j=0;j<NumberOfMolecules;j++) {
     345                        //if (TEList[j] != 0) {
     346                        runner = ListOfMolecules[j]->elemente->start;
     347                        while (runner->next != ListOfMolecules[j]->elemente->end) { // go through every element
     348                                runner = runner->next;
     349                                if (ListOfMolecules[j]->ElementsInMolecule[runner->Z]) { // if this element got atoms
     350                                        Walker = ListOfMolecules[j]->start;
     351                                        while (Walker->next != ListOfMolecules[j]->end) { // go through every atom of this element
     352                                                Walker = Walker->next;
     353                                                if (Walker->type->Z == runner->Z) {
     354                                                        if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
     355                                                                //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
     356                                                                ForcesFile <<   SortIndex[Walker->GetTrueFather()->nr] << "\t";
     357                                                                } else  // otherwise a -1 to indicate an added saturation hydrogen
     358                                                                        ForcesFile << "-1\t";
     359                                                        }
     360                                                }
     361                                        }
     362                        }
     363                        ForcesFile << endl;
     364                }
     365                ForcesFile.close();
     366                *out << Verbose(1) << "done." << endl;
     367        } else {
     368                status = false;
     369                *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
     370        }
     371        ForcesFile.close();
     372       
     373        return status;
    374374};
    375375
    … …  
    382382bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
    383383{
    384   ofstream outputFragment;
    385   char FragmentName[MAXSTRINGSIZE];
    386   char PathBackup[MAXSTRINGSIZE];
    387   bool result = true;
    388   bool intermediateResult = true;
    389   atom *Walker = NULL;
    390   Vector BoxDimension;
    391   char *FragmentNumber = NULL;
    392   char *path = NULL;
    393   int FragmentCounter = 0;
    394   ofstream output;
    395  
    396   // store the fragments as config and as xyz
    397   for(int i=0;i<NumberOfMolecules;i++) {
    398     // save default path as it is changed for each fragment
    399     path = configuration->GetDefaultPath();
    400     if (path != NULL)
    401       strcpy(PathBackup, path);
    402     else
    403       cerr << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl;
    404 
    405     // correct periodic
    406     ListOfMolecules[i]->ScanForPeriodicCorrection(out);
    407 
    408     // output xyz file
    409     FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
    410     sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    411     outputFragment.open(FragmentName, ios::out);
    412     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
    413     if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
    414       *out << " done." << endl;
    415     else
    416       *out << " failed." << endl;
    417     result = result && intermediateResult;
    418     outputFragment.close();
    419     outputFragment.clear();
    420 
    421     // list atoms in fragment for debugging
    422     *out << Verbose(2) << "Contained atoms: ";
    423     Walker = ListOfMolecules[i]->start;
    424     while (Walker->next != ListOfMolecules[i]->end) {
    425       Walker = Walker->next;
    426       *out << Walker->Name << " ";
    427     }
    428     *out << endl;
    429    
    430     // center on edge
    431     ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
    432     ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
    433     int j = -1;
    434     for (int k=0;k<NDIM;k++) {
    435       j += k+1;
    436       BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
    437       ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
    438     }
    439     ListOfMolecules[i]->Translate(&BoxDimension);
    440 
    441     // also calculate necessary orbitals
    442     ListOfMolecules[i]->CountElements();  // this is a bugfix, atoms should should actually be added correctly to this fragment
    443     ListOfMolecules[i]->CalculateOrbitals(*configuration);
    444    
    445     // change path in config
    446     //strcpy(PathBackup, configuration->configpath);
    447     sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
    448     configuration->SetDefaultPath(FragmentName);
    449    
    450     // and save as config
    451     sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    452     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
    453     if ((intermediateResult = configuration->Save(FragmentName, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
    454       *out << " done." << endl;
    455     else
    456       *out << " failed." << endl;
    457     result = result && intermediateResult;
    458 
    459     // restore old config
    460     configuration->SetDefaultPath(PathBackup);
    461 
    462 
    463     // and save as mpqc input file
    464     sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    465     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";
    466     if ((intermediateResult = configuration->SaveMPQC(FragmentName, ListOfMolecules[i])))
    467       *out << " done." << endl;
    468     else
    469       *out << " failed." << endl;
    470      
    471     result = result && intermediateResult;
    472     //outputFragment.close();
    473     //outputFragment.clear();
    474     delete(FragmentNumber);
    475     //Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
    476   }
    477   cout << " done." << endl;
    478  
    479   // printing final number
    480   *out << "Final number of fragments: " << FragmentCounter << "." << endl;
    481      
    482   return result;
     384        ofstream outputFragment;
     385        char FragmentName[MAXSTRINGSIZE];
     386        char PathBackup[MAXSTRINGSIZE];
     387        bool result = true;
     388        bool intermediateResult = true;
     389        atom *Walker = NULL;
     390        Vector BoxDimension;
     391        char *FragmentNumber = NULL;
     392        char *path = NULL;
     393        int FragmentCounter = 0;
     394        ofstream output;
     395       
     396        // store the fragments as config and as xyz
     397        for(int i=0;i<NumberOfMolecules;i++) {
     398                // save default path as it is changed for each fragment
     399                path = configuration->GetDefaultPath();
     400                if (path != NULL)
     401                        strcpy(PathBackup, path);
     402                else
     403                        cerr << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl;
     404
     405                // correct periodic
     406                ListOfMolecules[i]->ScanForPeriodicCorrection(out);
     407
     408                // output xyz file
     409                FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
     410                sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     411                outputFragment.open(FragmentName, ios::out);
     412                *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
     413                if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
     414                        *out << " done." << endl;
     415                else
     416                        *out << " failed." << endl;
     417                result = result && intermediateResult;
     418                outputFragment.close();
     419                outputFragment.clear();
     420
     421                // list atoms in fragment for debugging
     422                *out << Verbose(2) << "Contained atoms: ";
     423                Walker = ListOfMolecules[i]->start;
     424                while (Walker->next != ListOfMolecules[i]->end) {
     425                        Walker = Walker->next;
     426                        *out << Walker->Name << " ";
     427                }
     428                *out << endl;
     429               
     430                // center on edge
     431                ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
     432                ListOfMolecules[i]->SetBoxDimension(&BoxDimension);     // update Box of atoms by boundary
     433                int j = -1;
     434                for (int k=0;k<NDIM;k++) {
     435                        j += k+1;
     436                        BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
     437                        ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
     438                }
     439                ListOfMolecules[i]->Translate(&BoxDimension);
     440
     441                // also calculate necessary orbitals
     442                ListOfMolecules[i]->CountElements();    // this is a bugfix, atoms should should actually be added correctly to this fragment
     443                ListOfMolecules[i]->CalculateOrbitals(*configuration);
     444               
     445                // change path in config
     446                //strcpy(PathBackup, configuration->configpath);
     447                sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
     448                configuration->SetDefaultPath(FragmentName);
     449               
     450                // and save as config
     451                sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     452                *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
     453                if ((intermediateResult = configuration->Save(FragmentName, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
     454                        *out << " done." << endl;
     455                else
     456                        *out << " failed." << endl;
     457                result = result && intermediateResult;
     458
     459                // restore old config
     460                configuration->SetDefaultPath(PathBackup);
     461
     462
     463                // and save as mpqc input file
     464                sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     465                *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";
     466                if ((intermediateResult = configuration->SaveMPQC(FragmentName, ListOfMolecules[i])))
     467                        *out << " done." << endl;
     468                else
     469                        *out << " failed." << endl;
     470                       
     471                result = result && intermediateResult;
     472                //outputFragment.close();
     473                //outputFragment.clear();
     474                delete(FragmentNumber);
     475                //Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
     476        }
     477        cout << " done." << endl;
     478       
     479        // printing final number
     480        *out << "Final number of fragments: " << FragmentCounter << "." << endl;
     481                       
     482        return result;
    483483};
    484484
    … …  
    492492MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf = NULL)
    493493{
    494 //  if (Up != NULL)
    495 //    if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
    496 //      Up->DownLeaf = this;
    497 //  UpLeaf = Up;
    498 //  DownLeaf = NULL;
    499   Leaf = NULL;
    500   previous = PreviousLeaf;
    501   if (previous != NULL) {
    502     MoleculeLeafClass *Walker = previous->next;
    503     previous->next = this;
    504     next = Walker;
    505   } else {
    506     next = NULL;
    507   }
     494//      if (Up != NULL)
     495//              if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
     496//                      Up->DownLeaf = this;
     497//      UpLeaf = Up;
     498//      DownLeaf = NULL;
     499        Leaf = NULL;
     500        previous = PreviousLeaf;
     501        if (previous != NULL) {
     502                MoleculeLeafClass *Walker = previous->next;
     503                previous->next = this;
     504                next = Walker;
     505        } else {
     506                next = NULL;
     507        }
    508508};
    509509
    … …  
    512512MoleculeLeafClass::~MoleculeLeafClass()
    513513{
    514 //  if (DownLeaf != NULL) {// drop leaves further down
    515 //    MoleculeLeafClass *Walker = DownLeaf;
    516 //    MoleculeLeafClass *Next;
    517 //    do {
    518 //      Next = Walker->NextLeaf;
    519 //      delete(Walker);
    520 //      Walker = Next;
    521 //    } while (Walker != NULL);
    522 //    // Last Walker sets DownLeaf automatically to NULL
    523 //  }
    524   // remove the leaf itself
    525   if (Leaf != NULL) {
    526     delete(Leaf);
    527     Leaf = NULL;
    528   }
    529   // remove this Leaf from level list
    530   if (previous != NULL) 
    531     previous->next = next;
    532 //  } else { // we are first in list (connects to UpLeaf->DownLeaf)
    533 //    if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
    534 //      NextLeaf->UpLeaf = UpLeaf;  // either null as we are top level or the upleaf of the first node
    535 //    if (UpLeaf != NULL)
    536 //      UpLeaf->DownLeaf = NextLeaf;  // either null as we are only leaf or NextLeaf if we are just the first
    537 //  }
    538 //  UpLeaf = NULL;
    539   if (next != NULL) // are we last in list
    540     next->previous = previous;
    541   next = NULL;
    542   previous = NULL;
     514//      if (DownLeaf != NULL) {// drop leaves further down
     515//              MoleculeLeafClass *Walker = DownLeaf;
     516//              MoleculeLeafClass *Next;
     517//              do {
     518//                      Next = Walker->NextLeaf;
     519//                      delete(Walker);
     520//                      Walker = Next;
     521//              } while (Walker != NULL);
     522//              // Last Walker sets DownLeaf automatically to NULL
     523//      }
     524        // remove the leaf itself
     525        if (Leaf != NULL) {
     526                delete(Leaf);
     527                Leaf = NULL;
     528        }
     529        // remove this Leaf from level list
     530        if (previous != NULL)   
     531                previous->next = next;
     532//      } else { // we are first in list (connects to UpLeaf->DownLeaf)
     533//              if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
     534//                      NextLeaf->UpLeaf = UpLeaf;      // either null as we are top level or the upleaf of the first node
     535//              if (UpLeaf != NULL)
     536//                      UpLeaf->DownLeaf = NextLeaf;    // either null as we are only leaf or NextLeaf if we are just the first
     537//      }
     538//      UpLeaf = NULL;
     539        if (next != NULL) // are we last in list
     540                next->previous = previous;
     541        next = NULL;
     542        previous = NULL;
    543543};
    544544
    … …  
    550550bool MoleculeLeafClass::AddLeaf(molecule *ptr, MoleculeLeafClass *Previous)
    551551{
    552   return false;
     552        return false;
    553553};
    554554
    … …  
    564564bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
    565565{
    566   atom *Walker = NULL, *OtherWalker = NULL;
    567   bond *Binder = NULL;
    568   bool status = true;
    569   int AtomNo;
    570 
    571   *out << Verbose(1) << "Begin of FillBondStructureFromReference." << endl;
    572   // fill ListOfLocalAtoms if NULL was given
    573   if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
    574     *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
    575     return false;
    576   }
    577  
    578   if (status) {
    579     *out << Verbose(1) << "Creating adjacency list for subgraph " << this << "." << endl;
    580     Walker = Leaf->start;
    581     while (Walker->next != Leaf->end) {
    582       Walker = Walker->next;
    583       AtomNo = Walker->GetTrueFather()->nr;  // global id of the current walker
    584       for(int i=0;i<reference->NumberOfBondsPerAtom[AtomNo];i++) { // go through father's bonds and copy them all
    585         Binder = reference->ListOfBondsPerAtom[AtomNo][i];
    586         OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr];    // local copy of current bond partner of walker
    587         if (OtherWalker != NULL) {
    588           if (OtherWalker->nr > Walker->nr)
    589           Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
    590         } else {
    591           *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
    592           status = false;
    593         }
    594       }
    595     }
    596     Leaf->CreateListOfBondsPerAtom(out);
    597     FragmentCounter++;
    598     if (next != NULL)
    599       status = next->FillBondStructureFromReference(out, reference, FragmentCounter, ListOfLocalAtoms);
    600     FragmentCounter--;
    601   }
    602  
    603   if ((FreeList) && (ListOfLocalAtoms != NULL)) {
    604     // free the index lookup list
    605     Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
    606     if (FragmentCounter == 0) // first fragments frees the initial pointer to list
    607       Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
    608   }
    609   FragmentCounter--;
    610   *out << Verbose(1) << "End of FillBondStructureFromReference." << endl;
    611   return status;
     566        atom *Walker = NULL, *OtherWalker = NULL;
     567        bond *Binder = NULL;
     568        bool status = true;
     569        int AtomNo;
     570
     571        *out << Verbose(1) << "Begin of FillBondStructureFromReference." << endl;
     572        // fill ListOfLocalAtoms if NULL was given
     573        if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
     574                *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
     575                return false;
     576        }
     577       
     578        if (status) {
     579                *out << Verbose(1) << "Creating adjacency list for subgraph " << this << "." << endl;
     580                Walker = Leaf->start;
     581                while (Walker->next != Leaf->end) {
     582                        Walker = Walker->next;
     583                        AtomNo = Walker->GetTrueFather()->nr;   // global id of the current walker
     584                        for(int i=0;i<reference->NumberOfBondsPerAtom[AtomNo];i++) { // go through father's bonds and copy them all
     585                                Binder = reference->ListOfBondsPerAtom[AtomNo][i];
     586                                OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr];             // local copy of current bond partner of walker
     587                                if (OtherWalker != NULL) {
     588                                        if (OtherWalker->nr > Walker->nr)
     589                                        Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
     590                                } else {
     591                                        *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
     592                                        status = false;
     593                                }
     594                        }
     595                }
     596                Leaf->CreateListOfBondsPerAtom(out);
     597                FragmentCounter++;
     598                if (next != NULL)
     599                        status = next->FillBondStructureFromReference(out, reference, FragmentCounter, ListOfLocalAtoms);
     600                FragmentCounter--;
     601        }
     602       
     603        if ((FreeList) && (ListOfLocalAtoms != NULL)) {
     604                // free the index lookup list
     605                Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
     606                if (FragmentCounter == 0) // first fragments frees the initial pointer to list
     607                        Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
     608        }
     609        FragmentCounter--;
     610        *out << Verbose(1) << "End of FillBondStructureFromReference." << endl;
     611        return status;
    612612};
    613613
    … …  
    622622bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
    623623{
    624   atom *Walker = NULL, *Father = NULL;
    625 
    626   if (RootStack != NULL) {
    627     // find first root candidates
    628     if (&(RootStack[FragmentCounter]) != NULL) {
    629       RootStack[FragmentCounter].clear(); 
    630       Walker = Leaf->start;
    631       while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
    632         Walker = Walker->next;
    633         Father = Walker->GetTrueFather();
    634         if (AtomMask[Father->nr]) // apply mask
    635     #ifdef ADDHYDROGEN
    636           if (Walker->type->Z != 1) // skip hydrogen
    637     #endif
    638             RootStack[FragmentCounter].push_front(Walker->nr);
    639       }
    640       if (next != NULL)
    641         next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter);
    642     }  else {
    643       *out << Verbose(1) << "Rootstack[" << FragmentCounter  << "] is NULL." << endl;
    644       return false;
    645     }
    646     FragmentCounter--;
    647     return true;
    648   } else {
    649     *out << Verbose(1) << "Rootstack is NULL." << endl;
    650     return false;
    651   }
     624        atom *Walker = NULL, *Father = NULL;
     625
     626        if (RootStack != NULL) {
     627                // find first root candidates
     628                if (&(RootStack[FragmentCounter]) != NULL) {
     629                        RootStack[FragmentCounter].clear();     
     630                        Walker = Leaf->start;
     631                        while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
     632                                Walker = Walker->next;
     633                                Father = Walker->GetTrueFather();
     634                                if (AtomMask[Father->nr]) // apply mask
     635                #ifdef ADDHYDROGEN
     636                                        if (Walker->type->Z != 1) // skip hydrogen
     637                #endif
     638                                                RootStack[FragmentCounter].push_front(Walker->nr);
     639                        }
     640                        if (next != NULL)
     641                                next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter);
     642                }       else {
     643                        *out << Verbose(1) << "Rootstack[" << FragmentCounter   << "] is NULL." << endl;
     644                        return false;
     645                }
     646                FragmentCounter--;
     647                return true;
     648        } else {
     649                *out << Verbose(1) << "Rootstack is NULL." << endl;
     650                return false;
     651        }
    652652};
    653653
    … …  
    662662bool MoleculeLeafClass::FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList)
    663663{
    664   bool status = true;
    665  
    666   int Counter = Count();
    667   if (ListOfLocalAtoms == NULL) { // allocated initial pointer
    668     // allocate and set each field to NULL
    669     ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **)*Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
    670     if (ListOfLocalAtoms != NULL) {
    671       for (int i=Counter;i--;)
    672         ListOfLocalAtoms[i] = NULL;
    673       FreeList = FreeList && true;
    674     } else
    675       status = false;
    676   }
    677  
    678   if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
    679     status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
    680     FreeList = FreeList && true;
    681   }
    682  
    683   return status;
     664        bool status = true;
     665       
     666        int Counter = Count();
     667        if (ListOfLocalAtoms == NULL) { // allocated initial pointer
     668                // allocate and set each field to NULL
     669                ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **)*Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
     670                if (ListOfLocalAtoms != NULL) {
     671                        for (int i=Counter;i--;)
     672                                ListOfLocalAtoms[i] = NULL;
     673                        FreeList = FreeList && true;
     674                } else
     675                        status = false;
     676        }
     677       
     678        if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
     679                status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
     680                FreeList = FreeList && true;
     681        }
     682       
     683        return status;
    684684};
    685685
    … …  
    696696bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList)
    697697{
    698   bool status = true;
    699   int KeySetCounter = 0;
    700  
    701   *out << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl;
    702   // fill ListOfLocalAtoms if NULL was given
    703   if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
    704     *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
    705     return false;
    706   }
    707 
    708   // allocate fragment list
    709   if (FragmentList == NULL) {
    710     KeySetCounter = Count();
    711     FragmentList = (Graph **) Malloc(sizeof(Graph *)*KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
    712     for(int i=KeySetCounter;i--;)
    713       FragmentList[i] = NULL;
    714     KeySetCounter = 0;
    715   }
    716  
    717   if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
    718     // assign scanned keysets
    719     if (FragmentList[FragmentCounter] == NULL)
    720       FragmentList[FragmentCounter] = new Graph;
    721     KeySet *TempSet = new KeySet;
    722     for(Graph::iterator runner = KeySetList->begin();runner != KeySetList->end(); runner++) { // key sets contain global numbers!
    723       if ( ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
    724         // translate keyset to local numbers
    725         for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
    726           TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
    727         // insert into FragmentList
    728         FragmentList[FragmentCounter]->insert(GraphPair (*TempSet, pair<int,double>(KeySetCounter++, (*runner).second.second)));
    729       }
    730       TempSet->clear();
    731     }
    732     delete(TempSet);
    733     if (KeySetCounter == 0) {// if there are no keysets, delete the list
    734       *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
    735       delete(FragmentList[FragmentCounter]);
    736     } else
    737       *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
    738     FragmentCounter++;
    739     if (next != NULL)
    740       next->AssignKeySetsToFragment(out, reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
    741     FragmentCounter--;
    742   } else
    743     *out << Verbose(1) << "KeySetList is NULL or empty." << endl;
    744  
    745   if ((FreeList) && (ListOfLocalAtoms != NULL)) {
    746     // free the index lookup list
    747     Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]");
    748     if (FragmentCounter == 0) // first fragments frees the initial pointer to list
    749       Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms");
    750   }
    751   *out << Verbose(1) << "End of AssignKeySetsToFragment." << endl;
    752   return status;
     698        bool status = true;
     699        int KeySetCounter = 0;
     700       
     701        *out << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl;
     702        // fill ListOfLocalAtoms if NULL was given
     703        if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
     704                *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
     705                return false;
     706        }
     707
     708        // allocate fragment list
     709        if (FragmentList == NULL) {
     710                KeySetCounter = Count();
     711                FragmentList = (Graph **) Malloc(sizeof(Graph *)*KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
     712                for(int i=KeySetCounter;i--;)
     713                        FragmentList[i] = NULL;
     714                KeySetCounter = 0;
     715        }
     716       
     717        if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
     718                // assign scanned keysets
     719                if (FragmentList[FragmentCounter] == NULL)
     720                        FragmentList[FragmentCounter] = new Graph;
     721                KeySet *TempSet = new KeySet;
     722                for(Graph::iterator runner = KeySetList->begin();runner != KeySetList->end(); runner++) { // key sets contain global numbers!
     723                        if ( ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
     724                                // translate keyset to local numbers
     725                                for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
     726                                        TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
     727                                // insert into FragmentList
     728                                FragmentList[FragmentCounter]->insert(GraphPair (*TempSet, pair<int,double>(KeySetCounter++, (*runner).second.second)));
     729                        }
     730                        TempSet->clear();
     731                }
     732                delete(TempSet);
     733                if (KeySetCounter == 0) {// if there are no keysets, delete the list
     734                        *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
     735                        delete(FragmentList[FragmentCounter]);
     736                } else
     737                        *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
     738                FragmentCounter++;
     739                if (next != NULL)
     740                        next->AssignKeySetsToFragment(out, reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
     741                FragmentCounter--;
     742        } else
     743                *out << Verbose(1) << "KeySetList is NULL or empty." << endl;
     744       
     745        if ((FreeList) && (ListOfLocalAtoms != NULL)) {
     746                // free the index lookup list
     747                Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]");
     748                if (FragmentCounter == 0) // first fragments frees the initial pointer to list
     749                        Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms");
     750        }
     751        *out << Verbose(1) << "End of AssignKeySetsToFragment." << endl;
     752        return status;
    753753};
    754754
    … …  
    762762void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph)
    763763{
    764   *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
    765   KeySet *TempSet = new KeySet;
    766   if (FragmentList[FragmentCounter] != NULL) {
    767     for(Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
    768       for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
    769         TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
    770       TotalGraph.insert(GraphPair(*TempSet, pair<int,double>(TotalNumberOfKeySets++, (*runner).second.second)));
    771       TempSet->clear();
    772     }
    773     delete(TempSet);
    774   } else {
    775     *out << Verbose(1) << "FragmentList is NULL." << endl;
    776   }
    777   if (next != NULL)
    778     next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
    779   FragmentCounter--;
    780   *out << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl;
     764        *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
     765        KeySet *TempSet = new KeySet;
     766        if (FragmentList[FragmentCounter] != NULL) {
     767                for(Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
     768                        for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
     769                                TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
     770                        TotalGraph.insert(GraphPair(*TempSet, pair<int,double>(TotalNumberOfKeySets++, (*runner).second.second)));
     771                        TempSet->clear();
     772                }
     773                delete(TempSet);
     774        } else {
     775                *out << Verbose(1) << "FragmentList is NULL." << endl;
     776        }
     777        if (next != NULL)
     778                next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
     779        FragmentCounter--;
     780        *out << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl;
    781781};
    782782
    … …  
    786786int MoleculeLeafClass::Count() const
    787787{
    788   if (next != NULL)
    789     return next->Count()+1;
    790   else
    791     return 1;
    792 };
     788        if (next != NULL)
     789                return next->Count()+1;
     790        else
     791                return 1;
     792};
Note: See TracChangeset for help on using the changeset viewer.