1 | /** \file parsing.cpp
|
---|
2 | *
|
---|
3 | * Declarations of various class functions for the parsing of value files.
|
---|
4 | *
|
---|
5 | */
|
---|
6 |
|
---|
7 | // ======================================= INCLUDES ==========================================
|
---|
8 |
|
---|
9 | #include "Helpers/MemDebug.hpp"
|
---|
10 |
|
---|
11 | #include <cstring>
|
---|
12 |
|
---|
13 | #include "helpers.hpp"
|
---|
14 | #include "memoryallocator.hpp"
|
---|
15 | #include "parser.hpp"
|
---|
16 |
|
---|
17 | // include config.h
|
---|
18 | #ifdef HAVE_CONFIG_H
|
---|
19 | #include <config.h>
|
---|
20 | #endif
|
---|
21 |
|
---|
22 | // ======================================= FUNCTIONS ==========================================
|
---|
23 |
|
---|
24 | /** Test if given filename can be opened.
|
---|
25 | * \param filename name of file
|
---|
26 | * \param test true - no error message, false - print error
|
---|
27 | * \return given file exists
|
---|
28 | */
|
---|
29 | bool FilePresent(const char *filename, bool test)
|
---|
30 | {
|
---|
31 | ifstream input;
|
---|
32 |
|
---|
33 | input.open(filename, ios::in);
|
---|
34 | if (input == NULL) {
|
---|
35 | if (!test)
|
---|
36 | DoLog(0) && (Log() << Verbose(0) << endl << "FilePresent: Unable to open " << filename << ", is the directory correct?" << endl);
|
---|
37 | return false;
|
---|
38 | }
|
---|
39 | input.close();
|
---|
40 | return true;
|
---|
41 | };
|
---|
42 |
|
---|
43 | /** Test the given parameters.
|
---|
44 | * \param argc argument count
|
---|
45 | * \param **argv arguments array
|
---|
46 | * \return given inputdir is valid
|
---|
47 | */
|
---|
48 | bool TestParams(int argc, char **argv)
|
---|
49 | {
|
---|
50 | ifstream input;
|
---|
51 | stringstream line;
|
---|
52 |
|
---|
53 | line << argv[1] << FRAGMENTPREFIX << KEYSETFILE;
|
---|
54 | return FilePresent(line.str().c_str(), false);
|
---|
55 | };
|
---|
56 |
|
---|
57 | // ======================================= CLASS MatrixContainer =============================
|
---|
58 |
|
---|
59 | /** Constructor of MatrixContainer class.
|
---|
60 | */
|
---|
61 | MatrixContainer::MatrixContainer() {
|
---|
62 | Indices = NULL;
|
---|
63 | Header = new char*[1];
|
---|
64 | Matrix = new double**[1]; // one more each for the total molecule
|
---|
65 | RowCounter = new int[1];
|
---|
66 | ColumnCounter = new int[1];
|
---|
67 | Header[0] = NULL;
|
---|
68 | Matrix[0] = NULL;
|
---|
69 | RowCounter[0] = -1;
|
---|
70 | MatrixCounter = 0;
|
---|
71 | ColumnCounter[0] = -1;
|
---|
72 | };
|
---|
73 |
|
---|
74 | /** Destructor of MatrixContainer class.
|
---|
75 | */
|
---|
76 | MatrixContainer::~MatrixContainer() {
|
---|
77 | if (Matrix != NULL) {
|
---|
78 | // free Matrix[MatrixCounter]
|
---|
79 | if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
|
---|
80 | for(int j=RowCounter[MatrixCounter]+1;j--;)
|
---|
81 | delete[](Matrix[MatrixCounter][j]);
|
---|
82 | //if (MatrixCounter != 0)
|
---|
83 | delete[](Matrix[MatrixCounter]);
|
---|
84 | // free all matrices but ultimate one
|
---|
85 | for(int i=MatrixCounter;i--;) {
|
---|
86 | if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
|
---|
87 | for(int j=RowCounter[i]+1;j--;)
|
---|
88 | delete[](Matrix[i][j]);
|
---|
89 | delete[](Matrix[i]);
|
---|
90 | }
|
---|
91 | }
|
---|
92 | delete[](Matrix);
|
---|
93 | }
|
---|
94 | // free indices
|
---|
95 | if (Indices != NULL)
|
---|
96 | for(int i=MatrixCounter+1;i--;) {
|
---|
97 | delete[](Indices[i]);
|
---|
98 | }
|
---|
99 | delete[](Indices);
|
---|
100 |
|
---|
101 | // free header and counters
|
---|
102 | if (Header != NULL)
|
---|
103 | for(int i=MatrixCounter+1;i--;)
|
---|
104 | delete[](Header[i]);
|
---|
105 | delete[](Header);
|
---|
106 | delete[](RowCounter);
|
---|
107 | delete[](ColumnCounter);
|
---|
108 | };
|
---|
109 |
|
---|
110 | /** Either copies index matrix from another MatrixContainer or initialises with trivial mapping if NULL.
|
---|
111 | * This either copies the index matrix or just maps 1 to 1, 2 to 2 and so on for all fragments.
|
---|
112 | * \param *Matrix pointer to other MatrixContainer
|
---|
113 | * \return true - copy/initialisation sucessful, false - dimension false for copying
|
---|
114 | */
|
---|
115 | bool MatrixContainer::InitialiseIndices(class MatrixContainer *Matrix)
|
---|
116 | {
|
---|
117 | DoLog(0) && (Log() << Verbose(0) << "Initialising indices");
|
---|
118 | if (Matrix == NULL) {
|
---|
119 | DoLog(0) && (Log() << Verbose(0) << " with trivial mapping." << endl);
|
---|
120 | Indices = new int*[MatrixCounter + 1];
|
---|
121 | for(int i=MatrixCounter+1;i--;) {
|
---|
122 | Indices[i] = new int[RowCounter[i]];
|
---|
123 | for(int j=RowCounter[i];j--;)
|
---|
124 | Indices[i][j] = j;
|
---|
125 | }
|
---|
126 | } else {
|
---|
127 | DoLog(0) && (Log() << Verbose(0) << " from other MatrixContainer." << endl);
|
---|
128 | if (MatrixCounter != Matrix->MatrixCounter)
|
---|
129 | return false;
|
---|
130 | Indices = new int*[MatrixCounter + 1];
|
---|
131 | for(int i=MatrixCounter+1;i--;) {
|
---|
132 | if (RowCounter[i] != Matrix->RowCounter[i])
|
---|
133 | return false;
|
---|
134 | Indices[i] = new int[Matrix->RowCounter[i]];
|
---|
135 | for(int j=Matrix->RowCounter[i];j--;) {
|
---|
136 | Indices[i][j] = Matrix->Indices[i][j];
|
---|
137 | //Log() << Verbose(0) << Indices[i][j] << "\t";
|
---|
138 | }
|
---|
139 | //Log() << Verbose(0) << endl;
|
---|
140 | }
|
---|
141 | }
|
---|
142 | return true;
|
---|
143 | };
|
---|
144 |
|
---|
145 | /** Parsing a number of matrices.
|
---|
146 | * -# open the matrix file
|
---|
147 | * -# skip some lines (\a skiplines)
|
---|
148 | * -# scan header lines for number of columns
|
---|
149 | * -# scan lines for number of rows
|
---|
150 | * -# allocate matrix
|
---|
151 | * -# loop over found column and row counts and parse in each entry
|
---|
152 | * \param *name directory with files
|
---|
153 | * \param skiplines number of inital lines to skip
|
---|
154 | * \param skiplines number of inital columns to skip
|
---|
155 | * \param MatrixNr index number in Matrix array to parse into
|
---|
156 | * \return parsing successful
|
---|
157 | */
|
---|
158 | bool MatrixContainer::ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr)
|
---|
159 | {
|
---|
160 | ifstream input;
|
---|
161 | stringstream line;
|
---|
162 | string token;
|
---|
163 | char filename[1023];
|
---|
164 |
|
---|
165 | input.open(name, ios::in);
|
---|
166 | //Log() << Verbose(1) << "Opening " << name << " ... " << input << endl;
|
---|
167 | if (input == NULL) {
|
---|
168 | DoeLog(1) && (eLog()<< Verbose(1) << endl << "MatrixContainer::ParseMatrix: Unable to open " << name << ", is the directory correct?" << endl);
|
---|
169 | //performCriticalExit();
|
---|
170 | return false;
|
---|
171 | }
|
---|
172 |
|
---|
173 | // parse header
|
---|
174 | if (Header[MatrixNr] != NULL)
|
---|
175 | delete[] Header[MatrixNr];
|
---|
176 | Header[MatrixNr] = new char[1024];
|
---|
177 | for (int m=skiplines+1;m--;)
|
---|
178 | input.getline(Header[MatrixNr], 1023);
|
---|
179 |
|
---|
180 | // scan header for number of columns
|
---|
181 | line.str(Header[MatrixNr]);
|
---|
182 | for(int k=skipcolumns;k--;)
|
---|
183 | line >> Header[MatrixNr];
|
---|
184 | //Log() << Verbose(0) << line.str() << endl;
|
---|
185 | ColumnCounter[MatrixNr]=0;
|
---|
186 | while ( getline(line,token, '\t') ) {
|
---|
187 | if (token.length() > 0)
|
---|
188 | ColumnCounter[MatrixNr]++;
|
---|
189 | }
|
---|
190 | //Log() << Verbose(0) << line.str() << endl;
|
---|
191 | //Log() << Verbose(1) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
|
---|
192 | if (ColumnCounter[MatrixNr] == 0) {
|
---|
193 | DoeLog(0) && (eLog()<< Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl);
|
---|
194 | performCriticalExit();
|
---|
195 | }
|
---|
196 |
|
---|
197 | // scan rest for number of rows/lines
|
---|
198 | RowCounter[MatrixNr]=-1; // counts one line too much
|
---|
199 | while (!input.eof()) {
|
---|
200 | input.getline(filename, 1023);
|
---|
201 | //Log() << Verbose(0) << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
|
---|
202 | RowCounter[MatrixNr]++; // then line was not last MeanForce
|
---|
203 | if (strncmp(filename,"MeanForce", 9) == 0) {
|
---|
204 | break;
|
---|
205 | }
|
---|
206 | }
|
---|
207 | //Log() << Verbose(1) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
|
---|
208 | if (RowCounter[MatrixNr] == 0) {
|
---|
209 | DoeLog(0) && (eLog()<< Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl);
|
---|
210 | performCriticalExit();
|
---|
211 | }
|
---|
212 |
|
---|
213 | // allocate matrix if it's not zero dimension in one direction
|
---|
214 | if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
|
---|
215 | if (Matrix[MatrixNr] != NULL)
|
---|
216 | delete[] Matrix[MatrixNr];
|
---|
217 | Matrix[MatrixNr] = new double*[RowCounter[MatrixNr] + 1];
|
---|
218 | for(int j=0;j<RowCounter[MatrixNr]+1;j++)
|
---|
219 | Matrix[MatrixNr][j] = 0;
|
---|
220 |
|
---|
221 | // parse in each entry for this matrix
|
---|
222 | input.clear();
|
---|
223 | input.seekg(ios::beg);
|
---|
224 | for (int m=skiplines+1;m--;)
|
---|
225 | input.getline(Header[MatrixNr], 1023); // skip header
|
---|
226 | line.str(Header[MatrixNr]);
|
---|
227 | for(int k=skipcolumns;k--;) // skip columns in header too
|
---|
228 | line >> filename;
|
---|
229 | strncpy(Header[MatrixNr], line.str().c_str(), 1023);
|
---|
230 | for(int j=0;j<RowCounter[MatrixNr];j++) {
|
---|
231 | if (Matrix[MatrixNr][j] != NULL)
|
---|
232 | delete[] Matrix[MatrixNr][j];
|
---|
233 | Matrix[MatrixNr][j] = new double[ColumnCounter[MatrixNr]];
|
---|
234 | for(int k=0;k<ColumnCounter[MatrixNr];k++)
|
---|
235 | Matrix[MatrixNr][j][k] = 0;
|
---|
236 |
|
---|
237 | input.getline(filename, 1023);
|
---|
238 | stringstream lines(filename);
|
---|
239 | //Log() << Verbose(2) << "Matrix at level " << j << ":";// << filename << endl;
|
---|
240 | for(int k=skipcolumns;k--;)
|
---|
241 | lines >> filename;
|
---|
242 | for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
|
---|
243 | lines >> Matrix[MatrixNr][j][k];
|
---|
244 | //Log() << Verbose(1) << " " << setprecision(2) << Matrix[MatrixNr][j][k] << endl;
|
---|
245 | }
|
---|
246 | if (Matrix[MatrixNr][ RowCounter[MatrixNr] ] != NULL)
|
---|
247 | delete[] Matrix[MatrixNr][ RowCounter[MatrixNr] ];
|
---|
248 | Matrix[MatrixNr][ RowCounter[MatrixNr] ] = new double[ColumnCounter[MatrixNr]];
|
---|
249 | for(int j=ColumnCounter[MatrixNr];j--;)
|
---|
250 | Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
|
---|
251 | }
|
---|
252 | } else {
|
---|
253 | DoeLog(1) && (eLog()<< Verbose(1) << "Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl);
|
---|
254 | }
|
---|
255 | input.close();
|
---|
256 | return true;
|
---|
257 | };
|
---|
258 |
|
---|
259 | /** Parsing a number of matrices.
|
---|
260 | * -# First, count the number of matrices by counting lines in KEYSETFILE
|
---|
261 | * -# Then,
|
---|
262 | * -# construct the fragment number
|
---|
263 | * -# open the matrix file
|
---|
264 | * -# skip some lines (\a skiplines)
|
---|
265 | * -# scan header lines for number of columns
|
---|
266 | * -# scan lines for number of rows
|
---|
267 | * -# allocate matrix
|
---|
268 | * -# loop over found column and row counts and parse in each entry
|
---|
269 | * -# Finally, allocate one additional matrix (\a MatrixCounter) containing combined or temporary values
|
---|
270 | * \param *name directory with files
|
---|
271 | * \param *prefix prefix of each matrix file
|
---|
272 | * \param *suffix suffix of each matrix file
|
---|
273 | * \param skiplines number of inital lines to skip
|
---|
274 | * \param skiplines number of inital columns to skip
|
---|
275 | * \return parsing successful
|
---|
276 | */
|
---|
277 | bool MatrixContainer::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
|
---|
278 | {
|
---|
279 | char filename[1023];
|
---|
280 | ifstream input;
|
---|
281 | char *FragmentNumber = NULL;
|
---|
282 | stringstream file;
|
---|
283 | string token;
|
---|
284 |
|
---|
285 | // count the number of matrices
|
---|
286 | MatrixCounter = -1; // we count one too much
|
---|
287 | file << name << FRAGMENTPREFIX << KEYSETFILE;
|
---|
288 | input.open(file.str().c_str(), ios::in);
|
---|
289 | if (input == NULL) {
|
---|
290 | DoLog(0) && (Log() << Verbose(0) << endl << "MatrixContainer::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?" << endl);
|
---|
291 | return false;
|
---|
292 | }
|
---|
293 | while (!input.eof()) {
|
---|
294 | input.getline(filename, 1023);
|
---|
295 | stringstream zeile(filename);
|
---|
296 | MatrixCounter++;
|
---|
297 | }
|
---|
298 | input.close();
|
---|
299 | DoLog(0) && (Log() << Verbose(0) << "Determined " << MatrixCounter << " fragments." << endl);
|
---|
300 |
|
---|
301 | DoLog(0) && (Log() << Verbose(0) << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl);
|
---|
302 | delete[](Header);
|
---|
303 | delete[](Matrix);
|
---|
304 | delete[](RowCounter);
|
---|
305 | delete[](ColumnCounter);
|
---|
306 | Header = new char*[MatrixCounter + 1]; // one more each for the total molecule
|
---|
307 | Matrix = new double**[MatrixCounter + 1]; // one more each for the total molecule
|
---|
308 | RowCounter = new int[MatrixCounter + 1];
|
---|
309 | ColumnCounter = new int[MatrixCounter + 1];
|
---|
310 | for(int i=MatrixCounter+1;i--;) {
|
---|
311 | Matrix[i] = NULL;
|
---|
312 | Header[i] = NULL;
|
---|
313 | RowCounter[i] = -1;
|
---|
314 | ColumnCounter[i] = -1;
|
---|
315 | }
|
---|
316 | for(int i=0; i < MatrixCounter;i++) {
|
---|
317 | // open matrix file
|
---|
318 | FragmentNumber = FixedDigitNumber(MatrixCounter, i);
|
---|
319 | file.str(" ");
|
---|
320 | file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
|
---|
321 | if (!ParseMatrix(file.str().c_str(), skiplines, skipcolumns, i))
|
---|
322 | return false;
|
---|
323 | delete[](FragmentNumber);
|
---|
324 | }
|
---|
325 | return true;
|
---|
326 | };
|
---|
327 |
|
---|
328 | /** Allocates and resets the memory for a number \a MCounter of matrices.
|
---|
329 | * \param **GivenHeader Header line for each matrix
|
---|
330 | * \param MCounter number of matrices
|
---|
331 | * \param *RCounter number of rows for each matrix
|
---|
332 | * \param *CCounter number of columns for each matrix
|
---|
333 | * \return Allocation successful
|
---|
334 | */
|
---|
335 | bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
|
---|
336 | {
|
---|
337 | MatrixCounter = MCounter;
|
---|
338 | Header = new char*[MatrixCounter + 1];
|
---|
339 | Matrix = new double**[MatrixCounter + 1]; // one more each for the total molecule
|
---|
340 | RowCounter = new int[MatrixCounter + 1];
|
---|
341 | ColumnCounter = new int[MatrixCounter + 1];
|
---|
342 | for(int i=MatrixCounter+1;i--;) {
|
---|
343 | Header[i] = new char[1024];
|
---|
344 | strncpy(Header[i], GivenHeader[i], 1023);
|
---|
345 | RowCounter[i] = RCounter[i];
|
---|
346 | ColumnCounter[i] = CCounter[i];
|
---|
347 | Matrix[i] = new double*[RowCounter[i] + 1];
|
---|
348 | for(int j=RowCounter[i]+1;j--;) {
|
---|
349 | Matrix[i][j] = new double[ColumnCounter[i]];
|
---|
350 | for(int k=ColumnCounter[i];k--;)
|
---|
351 | Matrix[i][j][k] = 0.;
|
---|
352 | }
|
---|
353 | }
|
---|
354 | return true;
|
---|
355 | };
|
---|
356 |
|
---|
357 | /** Resets all values in MatrixContainer::Matrix.
|
---|
358 | * \return true if successful
|
---|
359 | */
|
---|
360 | bool MatrixContainer::ResetMatrix()
|
---|
361 | {
|
---|
362 | for(int i=MatrixCounter+1;i--;)
|
---|
363 | for(int j=RowCounter[i]+1;j--;)
|
---|
364 | for(int k=ColumnCounter[i];k--;)
|
---|
365 | Matrix[i][j][k] = 0.;
|
---|
366 | return true;
|
---|
367 | };
|
---|
368 |
|
---|
369 | /** Scans all elements of MatrixContainer::Matrix for greatest absolute value.
|
---|
370 | * \return greatest value of MatrixContainer::Matrix
|
---|
371 | */
|
---|
372 | double MatrixContainer::FindMaxValue()
|
---|
373 | {
|
---|
374 | double max = Matrix[0][0][0];
|
---|
375 | for(int i=MatrixCounter+1;i--;)
|
---|
376 | for(int j=RowCounter[i]+1;j--;)
|
---|
377 | for(int k=ColumnCounter[i];k--;)
|
---|
378 | if (fabs(Matrix[i][j][k]) > max)
|
---|
379 | max = fabs(Matrix[i][j][k]);
|
---|
380 | if (fabs(max) < MYEPSILON)
|
---|
381 | max += MYEPSILON;
|
---|
382 | return max;
|
---|
383 | };
|
---|
384 |
|
---|
385 | /** Scans all elements of MatrixContainer::Matrix for smallest absolute value.
|
---|
386 | * \return smallest value of MatrixContainer::Matrix
|
---|
387 | */
|
---|
388 | double MatrixContainer::FindMinValue()
|
---|
389 | {
|
---|
390 | double min = Matrix[0][0][0];
|
---|
391 | for(int i=MatrixCounter+1;i--;)
|
---|
392 | for(int j=RowCounter[i]+1;j--;)
|
---|
393 | for(int k=ColumnCounter[i];k--;)
|
---|
394 | if (fabs(Matrix[i][j][k]) < min)
|
---|
395 | min = fabs(Matrix[i][j][k]);
|
---|
396 | if (fabs(min) < MYEPSILON)
|
---|
397 | min += MYEPSILON;
|
---|
398 | return min;
|
---|
399 | };
|
---|
400 |
|
---|
401 | /** Sets all values in the last of MatrixContainer::Matrix to \a value.
|
---|
402 | * \param value reset value
|
---|
403 | * \param skipcolumns skip initial columns
|
---|
404 | * \return true if successful
|
---|
405 | */
|
---|
406 | bool MatrixContainer::SetLastMatrix(double value, int skipcolumns)
|
---|
407 | {
|
---|
408 | for(int j=RowCounter[MatrixCounter]+1;j--;)
|
---|
409 | for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
|
---|
410 | Matrix[MatrixCounter][j][k] = value;
|
---|
411 | return true;
|
---|
412 | };
|
---|
413 |
|
---|
414 | /** Sets all values in the last of MatrixContainer::Matrix to \a value.
|
---|
415 | * \param **values matrix with each value (must have at least same dimensions!)
|
---|
416 | * \param skipcolumns skip initial columns
|
---|
417 | * \return true if successful
|
---|
418 | */
|
---|
419 | bool MatrixContainer::SetLastMatrix(double **values, int skipcolumns)
|
---|
420 | {
|
---|
421 | for(int j=RowCounter[MatrixCounter]+1;j--;)
|
---|
422 | for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
|
---|
423 | Matrix[MatrixCounter][j][k] = values[j][k];
|
---|
424 | return true;
|
---|
425 | };
|
---|
426 |
|
---|
427 | /** Sums the entries with each factor and put into last element of \a ***Matrix.
|
---|
428 | * Sums over "E"-terms to create the "F"-terms
|
---|
429 | * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
|
---|
430 | * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
|
---|
431 | * \param Order bond order
|
---|
432 | * \return true if summing was successful
|
---|
433 | */
|
---|
434 | bool MatrixContainer::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order)
|
---|
435 | {
|
---|
436 | // go through each order
|
---|
437 | for (int CurrentFragment=0;CurrentFragment<KeySets.FragmentsPerOrder[Order];CurrentFragment++) {
|
---|
438 | //Log() << Verbose(0) << "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
|
---|
439 | // then go per order through each suborder and pick together all the terms that contain this fragment
|
---|
440 | for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
|
---|
441 | for (int j=0;j<KeySets.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
|
---|
442 | if (KeySets.Contains(KeySets.OrderSet[Order][CurrentFragment], KeySets.OrderSet[SubOrder][j])) {
|
---|
443 | //Log() << Verbose(0) << "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << "." << endl;
|
---|
444 | // if the fragment's indices are all in the current fragment
|
---|
445 | for(int k=0;k<RowCounter[ KeySets.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
|
---|
446 | int m = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][k];
|
---|
447 | //Log() << Verbose(0) << "Current index is " << k << "/" << m << "." << endl;
|
---|
448 | if (m != -1) { // if it's not an added hydrogen
|
---|
449 | for (int l=0;l<RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
|
---|
450 | //Log() << Verbose(0) << "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
|
---|
451 | if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
|
---|
452 | m = l;
|
---|
453 | break;
|
---|
454 | }
|
---|
455 | }
|
---|
456 | //Log() << Verbose(0) << "Corresponding index in CurrentFragment is " << m << "." << endl;
|
---|
457 | if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
|
---|
458 | DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment] << " current force index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
|
---|
459 | performCriticalExit();
|
---|
460 | return false;
|
---|
461 | }
|
---|
462 | if (Order == SubOrder) { // equal order is always copy from Energies
|
---|
463 | for(int l=ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column
|
---|
464 | Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
|
---|
465 | } else {
|
---|
466 | for(int l=ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l--;)
|
---|
467 | Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
|
---|
468 | }
|
---|
469 | }
|
---|
470 | //if ((ColumnCounter[ KeySets.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
|
---|
471 | //Log() << Verbose(0) << "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " << Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
|
---|
472 | }
|
---|
473 | } else {
|
---|
474 | //Log() << Verbose(0) << "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
|
---|
475 | }
|
---|
476 | }
|
---|
477 | }
|
---|
478 | //Log() << Verbose(0) << "Final Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << KeySets.AtomCounter[0]-1 << "][" << 1 << "] = " << Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][KeySets.AtomCounter[0]-1][1] << endl;
|
---|
479 | }
|
---|
480 |
|
---|
481 | return true;
|
---|
482 | };
|
---|
483 |
|
---|
484 | /** Writes the summed total fragment terms \f$F_{ij}\f$ to file.
|
---|
485 | * \param *name inputdir
|
---|
486 | * \param *prefix prefix before \a EnergySuffix
|
---|
487 | * \return file was written
|
---|
488 | */
|
---|
489 | bool MatrixContainer::WriteTotalFragments(const char *name, const char *prefix)
|
---|
490 | {
|
---|
491 | ofstream output;
|
---|
492 | char *FragmentNumber = NULL;
|
---|
493 |
|
---|
494 | DoLog(0) && (Log() << Verbose(0) << "Writing fragment files." << endl);
|
---|
495 | for(int i=0;i<MatrixCounter;i++) {
|
---|
496 | stringstream line;
|
---|
497 | FragmentNumber = FixedDigitNumber(MatrixCounter, i);
|
---|
498 | line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
|
---|
499 | delete[](FragmentNumber);
|
---|
500 | output.open(line.str().c_str(), ios::out);
|
---|
501 | if (output == NULL) {
|
---|
502 | DoeLog(0) && (eLog()<< Verbose(0) << "MatrixContainer::WriteTotalFragments: Unable to open output energy file " << line.str() << "!" << endl);
|
---|
503 | performCriticalExit();
|
---|
504 | return false;
|
---|
505 | }
|
---|
506 | output << Header[i] << endl;
|
---|
507 | for(int j=0;j<RowCounter[i];j++) {
|
---|
508 | for(int k=0;k<ColumnCounter[i];k++)
|
---|
509 | output << scientific << Matrix[i][j][k] << "\t";
|
---|
510 | output << endl;
|
---|
511 | }
|
---|
512 | output.close();
|
---|
513 | }
|
---|
514 | return true;
|
---|
515 | };
|
---|
516 |
|
---|
517 | /** Writes the summed total values in the last matrix to file.
|
---|
518 | * \param *name inputdir
|
---|
519 | * \param *prefix prefix
|
---|
520 | * \param *suffix suffix
|
---|
521 | * \return file was written
|
---|
522 | */
|
---|
523 | bool MatrixContainer::WriteLastMatrix(const char *name, const char *prefix, const char *suffix)
|
---|
524 | {
|
---|
525 | ofstream output;
|
---|
526 | stringstream line;
|
---|
527 |
|
---|
528 | DoLog(0) && (Log() << Verbose(0) << "Writing matrix values of " << suffix << "." << endl);
|
---|
529 | line << name << prefix << suffix;
|
---|
530 | output.open(line.str().c_str(), ios::out);
|
---|
531 | if (output == NULL) {
|
---|
532 | DoeLog(0) && (eLog()<< Verbose(0) << "MatrixContainer::WriteLastMatrix: Unable to open output matrix file " << line.str() << "!" << endl);
|
---|
533 | performCriticalExit();
|
---|
534 | return false;
|
---|
535 | }
|
---|
536 | output << Header[MatrixCounter] << endl;
|
---|
537 | for(int j=0;j<RowCounter[MatrixCounter];j++) {
|
---|
538 | for(int k=0;k<ColumnCounter[MatrixCounter];k++)
|
---|
539 | output << scientific << Matrix[MatrixCounter][j][k] << "\t";
|
---|
540 | output << endl;
|
---|
541 | }
|
---|
542 | output.close();
|
---|
543 | return true;
|
---|
544 | };
|
---|
545 |
|
---|
546 | // ======================================= CLASS EnergyMatrix =============================
|
---|
547 |
|
---|
548 | /** Create a trivial energy index mapping.
|
---|
549 | * This just maps 1 to 1, 2 to 2 and so on for all fragments.
|
---|
550 | * \return creation sucessful
|
---|
551 | */
|
---|
552 | bool EnergyMatrix::ParseIndices()
|
---|
553 | {
|
---|
554 | DoLog(0) && (Log() << Verbose(0) << "Parsing energy indices." << endl);
|
---|
555 | Indices = new int*[MatrixCounter + 1];
|
---|
556 | for(int i=MatrixCounter+1;i--;) {
|
---|
557 | Indices[i] = new int[RowCounter[i]];
|
---|
558 | for(int j=RowCounter[i];j--;)
|
---|
559 | Indices[i][j] = j;
|
---|
560 | }
|
---|
561 | return true;
|
---|
562 | };
|
---|
563 |
|
---|
564 | /** Sums the energy with each factor and put into last element of \a EnergyMatrix::Matrix.
|
---|
565 | * Sums over the "F"-terms in ANOVA decomposition.
|
---|
566 | * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
|
---|
567 | * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
|
---|
568 | * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
|
---|
569 | * \param Order bond order
|
---|
570 | * \parsm sign +1 or -1
|
---|
571 | * \return true if summing was successful
|
---|
572 | */
|
---|
573 | bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySets, int Order, double sign)
|
---|
574 | {
|
---|
575 | // sum energy
|
---|
576 | if (CorrectionFragments == NULL)
|
---|
577 | for(int i=KeySets.FragmentsPerOrder[Order];i--;)
|
---|
578 | for(int j=RowCounter[ KeySets.OrderSet[Order][i] ];j--;)
|
---|
579 | for(int k=ColumnCounter[ KeySets.OrderSet[Order][i] ];k--;)
|
---|
580 | Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySets.OrderSet[Order][i] ][j][k];
|
---|
581 | else
|
---|
582 | for(int i=KeySets.FragmentsPerOrder[Order];i--;)
|
---|
583 | for(int j=RowCounter[ KeySets.OrderSet[Order][i] ];j--;)
|
---|
584 | for(int k=ColumnCounter[ KeySets.OrderSet[Order][i] ];k--;)
|
---|
585 | Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySets.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySets.OrderSet[Order][i] ][j][k]);
|
---|
586 | return true;
|
---|
587 | };
|
---|
588 |
|
---|
589 | /** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
|
---|
590 | * \param *name directory with files
|
---|
591 | * \param *prefix prefix of each matrix file
|
---|
592 | * \param *suffix suffix of each matrix file
|
---|
593 | * \param skiplines number of inital lines to skip
|
---|
594 | * \param skiplines number of inital columns to skip
|
---|
595 | * \return parsing successful
|
---|
596 | */
|
---|
597 | bool EnergyMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
|
---|
598 | {
|
---|
599 | char filename[1024];
|
---|
600 | bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
|
---|
601 |
|
---|
602 | if (status) {
|
---|
603 | // count maximum of columns
|
---|
604 | RowCounter[MatrixCounter] = 0;
|
---|
605 | ColumnCounter[MatrixCounter] = 0;
|
---|
606 | for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
|
---|
607 | if (RowCounter[j] > RowCounter[MatrixCounter])
|
---|
608 | RowCounter[MatrixCounter] = RowCounter[j];
|
---|
609 | if (ColumnCounter[j] > ColumnCounter[MatrixCounter]) // take maximum of all for last matrix
|
---|
610 | ColumnCounter[MatrixCounter] = ColumnCounter[j];
|
---|
611 | }
|
---|
612 | // allocate last plus one matrix
|
---|
613 | DoLog(0) && (Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl);
|
---|
614 | Matrix[MatrixCounter] = new double*[RowCounter[MatrixCounter] + 1];
|
---|
615 | for(int j=0;j<=RowCounter[MatrixCounter];j++)
|
---|
616 | Matrix[MatrixCounter][j] = new double[ColumnCounter[MatrixCounter]];
|
---|
617 |
|
---|
618 | // try independently to parse global energysuffix file if present
|
---|
619 | strncpy(filename, name, 1023);
|
---|
620 | strncat(filename, prefix, 1023-strlen(filename));
|
---|
621 | strncat(filename, suffix.c_str(), 1023-strlen(filename));
|
---|
622 | ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
|
---|
623 | }
|
---|
624 | return status;
|
---|
625 | };
|
---|
626 |
|
---|
627 | // ======================================= CLASS ForceMatrix =============================
|
---|
628 |
|
---|
629 | /** Parsing force Indices of each fragment
|
---|
630 | * \param *name directory with \a ForcesFile
|
---|
631 | * \return parsing successful
|
---|
632 | */
|
---|
633 | bool ForceMatrix::ParseIndices(const char *name)
|
---|
634 | {
|
---|
635 | ifstream input;
|
---|
636 | char *FragmentNumber = NULL;
|
---|
637 | char filename[1023];
|
---|
638 | stringstream line;
|
---|
639 |
|
---|
640 | DoLog(0) && (Log() << Verbose(0) << "Parsing force indices for " << MatrixCounter << " matrices." << endl);
|
---|
641 | Indices = new int*[MatrixCounter + 1];
|
---|
642 | line << name << FRAGMENTPREFIX << FORCESFILE;
|
---|
643 | input.open(line.str().c_str(), ios::in);
|
---|
644 | //Log() << Verbose(0) << "Opening " << line.str() << " ... " << input << endl;
|
---|
645 | if (input == NULL) {
|
---|
646 | DoLog(0) && (Log() << Verbose(0) << endl << "ForceMatrix::ParseIndices: Unable to open " << line.str() << ", is the directory correct?" << endl);
|
---|
647 | return false;
|
---|
648 | }
|
---|
649 | for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
|
---|
650 | // get the number of atoms for this fragment
|
---|
651 | input.getline(filename, 1023);
|
---|
652 | line.str(filename);
|
---|
653 | // parse the values
|
---|
654 | Indices[i] = new int[RowCounter[i]];
|
---|
655 | FragmentNumber = FixedDigitNumber(MatrixCounter, i);
|
---|
656 | //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
|
---|
657 | delete[](FragmentNumber);
|
---|
658 | for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
|
---|
659 | line >> Indices[i][j];
|
---|
660 | //Log() << Verbose(0) << " " << Indices[i][j];
|
---|
661 | }
|
---|
662 | //Log() << Verbose(0) << endl;
|
---|
663 | }
|
---|
664 | Indices[MatrixCounter] = new int[RowCounter[MatrixCounter]];
|
---|
665 | for(int j=RowCounter[MatrixCounter];j--;) {
|
---|
666 | Indices[MatrixCounter][j] = j;
|
---|
667 | }
|
---|
668 | input.close();
|
---|
669 | return true;
|
---|
670 | };
|
---|
671 |
|
---|
672 |
|
---|
673 | /** Sums the forces and puts into last element of \a ForceMatrix::Matrix.
|
---|
674 | * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
|
---|
675 | * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
|
---|
676 | * \param Order bond order
|
---|
677 | * \param sign +1 or -1
|
---|
678 | * \return true if summing was successful
|
---|
679 | */
|
---|
680 | bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign)
|
---|
681 | {
|
---|
682 | int FragmentNr;
|
---|
683 | // sum forces
|
---|
684 | for(int i=0;i<KeySets.FragmentsPerOrder[Order];i++) {
|
---|
685 | FragmentNr = KeySets.OrderSet[Order][i];
|
---|
686 | for(int l=0;l<RowCounter[ FragmentNr ];l++) {
|
---|
687 | int j = Indices[ FragmentNr ][l];
|
---|
688 | if (j > RowCounter[MatrixCounter]) {
|
---|
689 | DoeLog(0) && (eLog()<< Verbose(0) << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl);
|
---|
690 | performCriticalExit();
|
---|
691 | return false;
|
---|
692 | }
|
---|
693 | if (j != -1) {
|
---|
694 | //if (j == 0) Log() << Verbose(0) << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
|
---|
695 | for(int k=2;k<ColumnCounter[MatrixCounter];k++)
|
---|
696 | Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
|
---|
697 | }
|
---|
698 | }
|
---|
699 | }
|
---|
700 | return true;
|
---|
701 | };
|
---|
702 |
|
---|
703 |
|
---|
704 | /** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
|
---|
705 | * \param *name directory with files
|
---|
706 | * \param *prefix prefix of each matrix file
|
---|
707 | * \param *suffix suffix of each matrix file
|
---|
708 | * \param skiplines number of inital lines to skip
|
---|
709 | * \param skiplines number of inital columns to skip
|
---|
710 | * \return parsing successful
|
---|
711 | */
|
---|
712 | bool ForceMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
|
---|
713 | {
|
---|
714 | char filename[1023];
|
---|
715 | ifstream input;
|
---|
716 | stringstream file;
|
---|
717 | int nr;
|
---|
718 | bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
|
---|
719 |
|
---|
720 | if (status) {
|
---|
721 | // count number of atoms for last plus one matrix
|
---|
722 | file << name << FRAGMENTPREFIX << KEYSETFILE;
|
---|
723 | input.open(file.str().c_str(), ios::in);
|
---|
724 | if (input == NULL) {
|
---|
725 | DoLog(0) && (Log() << Verbose(0) << endl << "ForceMatrix::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?" << endl);
|
---|
726 | return false;
|
---|
727 | }
|
---|
728 | RowCounter[MatrixCounter] = 0;
|
---|
729 | while (!input.eof()) {
|
---|
730 | input.getline(filename, 1023);
|
---|
731 | stringstream zeile(filename);
|
---|
732 | while (!zeile.eof()) {
|
---|
733 | zeile >> nr;
|
---|
734 | //Log() << Verbose(0) << "Current index: " << nr << "." << endl;
|
---|
735 | if (nr > RowCounter[MatrixCounter])
|
---|
736 | RowCounter[MatrixCounter] = nr;
|
---|
737 | }
|
---|
738 | }
|
---|
739 | RowCounter[MatrixCounter]++; // nr start at 0, count starts at 1
|
---|
740 | input.close();
|
---|
741 |
|
---|
742 | ColumnCounter[MatrixCounter] = 0;
|
---|
743 | for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
|
---|
744 | if (ColumnCounter[j] > ColumnCounter[MatrixCounter]) // take maximum of all for last matrix
|
---|
745 | ColumnCounter[MatrixCounter] = ColumnCounter[j];
|
---|
746 | }
|
---|
747 |
|
---|
748 | // allocate last plus one matrix
|
---|
749 | DoLog(0) && (Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl);
|
---|
750 | Matrix[MatrixCounter] = new double*[RowCounter[MatrixCounter] + 1];
|
---|
751 | for(int j=0;j<=RowCounter[MatrixCounter];j++)
|
---|
752 | Matrix[MatrixCounter][j] = new double[ColumnCounter[MatrixCounter]];
|
---|
753 |
|
---|
754 | // try independently to parse global forcesuffix file if present
|
---|
755 | strncpy(filename, name, 1023);
|
---|
756 | strncat(filename, prefix, 1023-strlen(filename));
|
---|
757 | strncat(filename, suffix.c_str(), 1023-strlen(filename));
|
---|
758 | ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
|
---|
759 | }
|
---|
760 |
|
---|
761 |
|
---|
762 | return status;
|
---|
763 | };
|
---|
764 |
|
---|
765 | // ======================================= CLASS HessianMatrix =============================
|
---|
766 |
|
---|
767 | /** Parsing force Indices of each fragment
|
---|
768 | * \param *name directory with \a ForcesFile
|
---|
769 | * \return parsing successful
|
---|
770 | */
|
---|
771 | bool HessianMatrix::ParseIndices(char *name)
|
---|
772 | {
|
---|
773 | ifstream input;
|
---|
774 | char *FragmentNumber = NULL;
|
---|
775 | char filename[1023];
|
---|
776 | stringstream line;
|
---|
777 |
|
---|
778 | DoLog(0) && (Log() << Verbose(0) << "Parsing hessian indices for " << MatrixCounter << " matrices." << endl);
|
---|
779 | Indices = new int*[MatrixCounter + 1];
|
---|
780 | line << name << FRAGMENTPREFIX << FORCESFILE;
|
---|
781 | input.open(line.str().c_str(), ios::in);
|
---|
782 | //Log() << Verbose(0) << "Opening " << line.str() << " ... " << input << endl;
|
---|
783 | if (input == NULL) {
|
---|
784 | DoLog(0) && (Log() << Verbose(0) << endl << "HessianMatrix::ParseIndices: Unable to open " << line.str() << ", is the directory correct?" << endl);
|
---|
785 | return false;
|
---|
786 | }
|
---|
787 | for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
|
---|
788 | // get the number of atoms for this fragment
|
---|
789 | input.getline(filename, 1023);
|
---|
790 | line.str(filename);
|
---|
791 | // parse the values
|
---|
792 | Indices[i] = new int[RowCounter[i]];
|
---|
793 | FragmentNumber = FixedDigitNumber(MatrixCounter, i);
|
---|
794 | //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
|
---|
795 | delete[](FragmentNumber);
|
---|
796 | for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
|
---|
797 | line >> Indices[i][j];
|
---|
798 | //Log() << Verbose(0) << " " << Indices[i][j];
|
---|
799 | }
|
---|
800 | //Log() << Verbose(0) << endl;
|
---|
801 | }
|
---|
802 | Indices[MatrixCounter] = new int[RowCounter[MatrixCounter]];
|
---|
803 | for(int j=RowCounter[MatrixCounter];j--;) {
|
---|
804 | Indices[MatrixCounter][j] = j;
|
---|
805 | }
|
---|
806 | input.close();
|
---|
807 | return true;
|
---|
808 | };
|
---|
809 |
|
---|
810 |
|
---|
811 | /** Sums the hessian entries and puts into last element of \a HessianMatrix::Matrix.
|
---|
812 | * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
|
---|
813 | * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
|
---|
814 | * \param Order bond order
|
---|
815 | * \param sign +1 or -1
|
---|
816 | * \return true if summing was successful
|
---|
817 | */
|
---|
818 | bool HessianMatrix::SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign)
|
---|
819 | {
|
---|
820 | int FragmentNr;
|
---|
821 | // sum forces
|
---|
822 | for(int i=0;i<KeySets.FragmentsPerOrder[Order];i++) {
|
---|
823 | FragmentNr = KeySets.OrderSet[Order][i];
|
---|
824 | for(int l=0;l<RowCounter[ FragmentNr ];l++) {
|
---|
825 | int j = Indices[ FragmentNr ][l];
|
---|
826 | if (j > RowCounter[MatrixCounter]) {
|
---|
827 | DoeLog(0) && (eLog()<< Verbose(0) << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl);
|
---|
828 | performCriticalExit();
|
---|
829 | return false;
|
---|
830 | }
|
---|
831 | if (j != -1) {
|
---|
832 | for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
|
---|
833 | int k = Indices[ FragmentNr ][m];
|
---|
834 | if (k > ColumnCounter[MatrixCounter]) {
|
---|
835 | DoeLog(0) && (eLog()<< Verbose(0) << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl);
|
---|
836 | performCriticalExit();
|
---|
837 | return false;
|
---|
838 | }
|
---|
839 | if (k != -1) {
|
---|
840 | //Log() << Verbose(0) << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
|
---|
841 | Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
|
---|
842 | }
|
---|
843 | }
|
---|
844 | }
|
---|
845 | }
|
---|
846 | }
|
---|
847 | return true;
|
---|
848 | };
|
---|
849 |
|
---|
850 | /** Constructor for class HessianMatrix.
|
---|
851 | */
|
---|
852 | HessianMatrix::HessianMatrix() : MatrixContainer()
|
---|
853 | {
|
---|
854 | IsSymmetric = true;
|
---|
855 | }
|
---|
856 |
|
---|
857 | /** Sums the hessian entries with each factor and put into last element of \a ***Matrix.
|
---|
858 | * Sums over "E"-terms to create the "F"-terms
|
---|
859 | * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
|
---|
860 | * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
|
---|
861 | * \param Order bond order
|
---|
862 | * \return true if summing was successful
|
---|
863 | */
|
---|
864 | bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order)
|
---|
865 | {
|
---|
866 | // go through each order
|
---|
867 | for (int CurrentFragment=0;CurrentFragment<KeySets.FragmentsPerOrder[Order];CurrentFragment++) {
|
---|
868 | //Log() << Verbose(0) << "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
|
---|
869 | // then go per order through each suborder and pick together all the terms that contain this fragment
|
---|
870 | for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
|
---|
871 | for (int j=0;j<KeySets.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
|
---|
872 | if (KeySets.Contains(KeySets.OrderSet[Order][CurrentFragment], KeySets.OrderSet[SubOrder][j])) {
|
---|
873 | //Log() << Verbose(0) << "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << "." << endl;
|
---|
874 | // if the fragment's indices are all in the current fragment
|
---|
875 | for(int k=0;k<RowCounter[ KeySets.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
|
---|
876 | int m = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][k];
|
---|
877 | //Log() << Verbose(0) << "Current row index is " << k << "/" << m << "." << endl;
|
---|
878 | if (m != -1) { // if it's not an added hydrogen
|
---|
879 | for (int l=0;l<RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
|
---|
880 | //Log() << Verbose(0) << "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
|
---|
881 | if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
|
---|
882 | m = l;
|
---|
883 | break;
|
---|
884 | }
|
---|
885 | }
|
---|
886 | //Log() << Verbose(0) << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
|
---|
887 | if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
|
---|
888 | DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment] << " current row index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
|
---|
889 | performCriticalExit();
|
---|
890 | return false;
|
---|
891 | }
|
---|
892 |
|
---|
893 | for(int l=0;l<ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l++) {
|
---|
894 | int n = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][l];
|
---|
895 | //Log() << Verbose(0) << "Current column index is " << l << "/" << n << "." << endl;
|
---|
896 | if (n != -1) { // if it's not an added hydrogen
|
---|
897 | for (int p=0;p<ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
|
---|
898 | //Log() << Verbose(0) << "Comparing " << n << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
|
---|
899 | if (n == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p]) {
|
---|
900 | n = p;
|
---|
901 | break;
|
---|
902 | }
|
---|
903 | }
|
---|
904 | //Log() << Verbose(0) << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
|
---|
905 | if (n > ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
|
---|
906 | DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment] << " current column index " << n << " is greater than " << ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
|
---|
907 | performCriticalExit();
|
---|
908 | return false;
|
---|
909 | }
|
---|
910 | if (Order == SubOrder) { // equal order is always copy from Energies
|
---|
911 | //Log() << Verbose(0) << "Adding " << MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
|
---|
912 | Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
|
---|
913 | } else {
|
---|
914 | //Log() << Verbose(0) << "Subtracting " << Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
|
---|
915 | Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
|
---|
916 | }
|
---|
917 | }
|
---|
918 | }
|
---|
919 | }
|
---|
920 | //if ((ColumnCounter[ KeySets.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
|
---|
921 | //Log() << Verbose(0) << "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " << Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
|
---|
922 | }
|
---|
923 | } else {
|
---|
924 | //Log() << Verbose(0) << "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
|
---|
925 | }
|
---|
926 | }
|
---|
927 | }
|
---|
928 | //Log() << Verbose(0) << "Final Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << KeySets.AtomCounter[0]-1 << "][" << 1 << "] = " << Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][KeySets.AtomCounter[0]-1][1] << endl;
|
---|
929 | }
|
---|
930 |
|
---|
931 | return true;
|
---|
932 | };
|
---|
933 |
|
---|
934 | /** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
|
---|
935 | * \param *name directory with files
|
---|
936 | * \param *prefix prefix of each matrix file
|
---|
937 | * \param *suffix suffix of each matrix file
|
---|
938 | * \param skiplines number of inital lines to skip
|
---|
939 | * \param skiplines number of inital columns to skip
|
---|
940 | * \return parsing successful
|
---|
941 | */
|
---|
942 | bool HessianMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
|
---|
943 | {
|
---|
944 | char filename[1023];
|
---|
945 | ifstream input;
|
---|
946 | stringstream file;
|
---|
947 | int nr;
|
---|
948 | bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
|
---|
949 |
|
---|
950 | if (status) {
|
---|
951 | // count number of atoms for last plus one matrix
|
---|
952 | file << name << FRAGMENTPREFIX << KEYSETFILE;
|
---|
953 | input.open(file.str().c_str(), ios::in);
|
---|
954 | if (input == NULL) {
|
---|
955 | DoLog(0) && (Log() << Verbose(0) << endl << "HessianMatrix::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?" << endl);
|
---|
956 | return false;
|
---|
957 | }
|
---|
958 | RowCounter[MatrixCounter] = 0;
|
---|
959 | ColumnCounter[MatrixCounter] = 0;
|
---|
960 | while (!input.eof()) {
|
---|
961 | input.getline(filename, 1023);
|
---|
962 | stringstream zeile(filename);
|
---|
963 | while (!zeile.eof()) {
|
---|
964 | zeile >> nr;
|
---|
965 | //Log() << Verbose(0) << "Current index: " << nr << "." << endl;
|
---|
966 | if (nr > RowCounter[MatrixCounter]) {
|
---|
967 | RowCounter[MatrixCounter] = nr;
|
---|
968 | ColumnCounter[MatrixCounter] = nr;
|
---|
969 | }
|
---|
970 | }
|
---|
971 | }
|
---|
972 | RowCounter[MatrixCounter]++; // nr start at 0, count starts at 1
|
---|
973 | ColumnCounter[MatrixCounter]++; // nr start at 0, count starts at 1
|
---|
974 | input.close();
|
---|
975 |
|
---|
976 | // allocate last plus one matrix
|
---|
977 | DoLog(0) && (Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl);
|
---|
978 | Matrix[MatrixCounter] = new double*[RowCounter[MatrixCounter] + 1];
|
---|
979 | for(int j=0;j<=RowCounter[MatrixCounter];j++)
|
---|
980 | Matrix[MatrixCounter][j] = new double[ColumnCounter[MatrixCounter]];
|
---|
981 |
|
---|
982 | // try independently to parse global forcesuffix file if present
|
---|
983 | strncpy(filename, name, 1023);
|
---|
984 | strncat(filename, prefix, 1023-strlen(filename));
|
---|
985 | strncat(filename, suffix.c_str(), 1023-strlen(filename));
|
---|
986 | ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
|
---|
987 | }
|
---|
988 |
|
---|
989 |
|
---|
990 | return status;
|
---|
991 | };
|
---|
992 |
|
---|
993 | // ======================================= CLASS KeySetsContainer =============================
|
---|
994 |
|
---|
995 | /** Constructor of KeySetsContainer class.
|
---|
996 | */
|
---|
997 | KeySetsContainer::KeySetsContainer() {
|
---|
998 | KeySets = NULL;
|
---|
999 | AtomCounter = NULL;
|
---|
1000 | FragmentCounter = 0;
|
---|
1001 | Order = 0;
|
---|
1002 | FragmentsPerOrder = 0;
|
---|
1003 | OrderSet = NULL;
|
---|
1004 | };
|
---|
1005 |
|
---|
1006 | /** Destructor of KeySetsContainer class.
|
---|
1007 | */
|
---|
1008 | KeySetsContainer::~KeySetsContainer() {
|
---|
1009 | for(int i=FragmentCounter;i--;)
|
---|
1010 | delete[](KeySets[i]);
|
---|
1011 | for(int i=Order;i--;)
|
---|
1012 | delete[](OrderSet[i]);
|
---|
1013 | delete[](KeySets);
|
---|
1014 | delete[](OrderSet);
|
---|
1015 | delete[](AtomCounter);
|
---|
1016 | delete[](FragmentsPerOrder);
|
---|
1017 | };
|
---|
1018 |
|
---|
1019 | /** Parsing KeySets into array.
|
---|
1020 | * \param *name directory with keyset file
|
---|
1021 | * \param *ACounter number of atoms per fragment
|
---|
1022 | * \param FCounter number of fragments
|
---|
1023 | * \return parsing succesful
|
---|
1024 | */
|
---|
1025 | bool KeySetsContainer::ParseKeySets(const char *name, const int *ACounter, const int FCounter) {
|
---|
1026 | ifstream input;
|
---|
1027 | char *FragmentNumber = NULL;
|
---|
1028 | stringstream file;
|
---|
1029 | char filename[1023];
|
---|
1030 |
|
---|
1031 | FragmentCounter = FCounter;
|
---|
1032 | DoLog(0) && (Log() << Verbose(0) << "Parsing key sets." << endl);
|
---|
1033 | KeySets = new int*[FragmentCounter];
|
---|
1034 | for(int i=FragmentCounter;i--;)
|
---|
1035 | KeySets[i] = NULL;
|
---|
1036 | file << name << FRAGMENTPREFIX << KEYSETFILE;
|
---|
1037 | input.open(file.str().c_str(), ios::in);
|
---|
1038 | if (input == NULL) {
|
---|
1039 | DoLog(0) && (Log() << Verbose(0) << endl << "KeySetsContainer::ParseKeySets: Unable to open " << file.str() << ", is the directory correct?" << endl);
|
---|
1040 | return false;
|
---|
1041 | }
|
---|
1042 |
|
---|
1043 | AtomCounter = new int[FragmentCounter];
|
---|
1044 | for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
|
---|
1045 | stringstream line;
|
---|
1046 | AtomCounter[i] = ACounter[i];
|
---|
1047 | // parse the values
|
---|
1048 | KeySets[i] = new int[AtomCounter[i]];
|
---|
1049 | for(int j=AtomCounter[i];j--;)
|
---|
1050 | KeySets[i][j] = -1;
|
---|
1051 | FragmentNumber = FixedDigitNumber(FragmentCounter, i);
|
---|
1052 | //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
|
---|
1053 | delete[](FragmentNumber);
|
---|
1054 | input.getline(filename, 1023);
|
---|
1055 | line.str(filename);
|
---|
1056 | for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
|
---|
1057 | line >> KeySets[i][j];
|
---|
1058 | //Log() << Verbose(0) << " " << KeySets[i][j];
|
---|
1059 | }
|
---|
1060 | //Log() << Verbose(0) << endl;
|
---|
1061 | }
|
---|
1062 | input.close();
|
---|
1063 | return true;
|
---|
1064 | };
|
---|
1065 |
|
---|
1066 | /** Parse many body terms, associating each fragment to a certain bond order.
|
---|
1067 | * \return parsing succesful
|
---|
1068 | */
|
---|
1069 | bool KeySetsContainer::ParseManyBodyTerms()
|
---|
1070 | {
|
---|
1071 | int Counter;
|
---|
1072 |
|
---|
1073 | DoLog(0) && (Log() << Verbose(0) << "Creating Fragment terms." << endl);
|
---|
1074 | // scan through all to determine maximum order
|
---|
1075 | Order=0;
|
---|
1076 | for(int i=FragmentCounter;i--;) {
|
---|
1077 | Counter=0;
|
---|
1078 | for(int j=AtomCounter[i];j--;)
|
---|
1079 | if (KeySets[i][j] != -1)
|
---|
1080 | Counter++;
|
---|
1081 | if (Counter > Order)
|
---|
1082 | Order = Counter;
|
---|
1083 | }
|
---|
1084 | DoLog(0) && (Log() << Verbose(0) << "Found Order is " << Order << "." << endl);
|
---|
1085 |
|
---|
1086 | // scan through all to determine fragments per order
|
---|
1087 | FragmentsPerOrder = new int[Order];
|
---|
1088 | for(int i=Order;i--;)
|
---|
1089 | FragmentsPerOrder[i] = 0;
|
---|
1090 | for(int i=FragmentCounter;i--;) {
|
---|
1091 | Counter=0;
|
---|
1092 | for(int j=AtomCounter[i];j--;)
|
---|
1093 | if (KeySets[i][j] != -1)
|
---|
1094 | Counter++;
|
---|
1095 | FragmentsPerOrder[Counter-1]++;
|
---|
1096 | }
|
---|
1097 | for(int i=0;i<Order;i++)
|
---|
1098 | DoLog(0) && (Log() << Verbose(0) << "Found No. of Fragments of Order " << i+1 << " is " << FragmentsPerOrder[i] << "." << endl);
|
---|
1099 |
|
---|
1100 | // scan through all to gather indices to each order set
|
---|
1101 | OrderSet = new int*[Order];
|
---|
1102 | for(int i=Order;i--;)
|
---|
1103 | OrderSet[i] = new int[FragmentsPerOrder[i]];
|
---|
1104 | for(int i=Order;i--;)
|
---|
1105 | FragmentsPerOrder[i] = 0;
|
---|
1106 | for(int i=FragmentCounter;i--;) {
|
---|
1107 | Counter=0;
|
---|
1108 | for(int j=AtomCounter[i];j--;)
|
---|
1109 | if (KeySets[i][j] != -1)
|
---|
1110 | Counter++;
|
---|
1111 | OrderSet[Counter-1][FragmentsPerOrder[Counter-1]] = i;
|
---|
1112 | FragmentsPerOrder[Counter-1]++;
|
---|
1113 | }
|
---|
1114 | DoLog(0) && (Log() << Verbose(0) << "Printing OrderSet." << endl);
|
---|
1115 | for(int i=0;i<Order;i++) {
|
---|
1116 | for (int j=0;j<FragmentsPerOrder[i];j++) {
|
---|
1117 | DoLog(0) && (Log() << Verbose(0) << " " << OrderSet[i][j]);
|
---|
1118 | }
|
---|
1119 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
1120 | }
|
---|
1121 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
1122 |
|
---|
1123 |
|
---|
1124 | return true;
|
---|
1125 | };
|
---|
1126 |
|
---|
1127 | /** Compares each entry in \a *SmallerSet if it is containted in \a *GreaterSet.
|
---|
1128 | * \param GreaterSet index to greater set
|
---|
1129 | * \param SmallerSet index to smaller set
|
---|
1130 | * \return true if all keys of SmallerSet contained in GreaterSet
|
---|
1131 | */
|
---|
1132 | bool KeySetsContainer::Contains(const int GreaterSet, const int SmallerSet)
|
---|
1133 | {
|
---|
1134 | bool result = true;
|
---|
1135 | bool intermediate;
|
---|
1136 | if ((GreaterSet < 0) || (SmallerSet < 0) || (GreaterSet > FragmentCounter) || (SmallerSet > FragmentCounter)) // index out of bounds
|
---|
1137 | return false;
|
---|
1138 | for(int i=AtomCounter[SmallerSet];i--;) {
|
---|
1139 | intermediate = false;
|
---|
1140 | for (int j=AtomCounter[GreaterSet];j--;)
|
---|
1141 | intermediate = (intermediate || ((KeySets[SmallerSet][i] == KeySets[GreaterSet][j]) || (KeySets[SmallerSet][i] == -1)));
|
---|
1142 | result = result && intermediate;
|
---|
1143 | }
|
---|
1144 |
|
---|
1145 | return result;
|
---|
1146 | };
|
---|
1147 |
|
---|
1148 |
|
---|
1149 | // ======================================= END =============================================
|
---|