source: src/parser.cpp@ bd6579

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since bd6579 was 3b15bb, checked in by Frederik Heber <heber@…>, 17 years ago

~MatrixContainer(): Bugfix - Free if Matrix[MatrixCounter] would also occur if nothing had been allocated so far

  • Property mode set to 100644
File size: 23.1 KB
Line 
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.hpp"
10#include "parser.hpp"
11
12// ======================================= FUNCTIONS ==========================================
13
14/** Test if given filename can be opened.
15 * \param filename name of file
16 * \return given file exists
17 */
18#ifdef HAVE_INLINE
19inline bool FilePresent(const char *filename)
20#else
21bool FilePresent(const char *filename)
22#endif
23{
24 ifstream input;
25
26 input.open(filename, ios::in);
27 if (input == NULL) {
28 cout << endl << "Unable to open " << filename << ", is the directory correct?" << endl;
29 return false;
30 }
31 input.close();
32 return true;
33};
34
35/** Test the given parameters.
36 * \param argc argument count
37 * \param **argv arguments array
38 * \return given inputdir is valid
39 */
40bool TestParams(int argc, char **argv)
41{
42 ifstream input;
43 stringstream line;
44
45 //line << argv[1] << FRAGMENTPREFIX << TEFACTORSFILE;
46 return FilePresent(line.str().c_str());
47};
48
49// ======================================= CLASS MatrixContainer =============================
50
51/** Constructor of MatrixContainer class.
52 */
53MatrixContainer::MatrixContainer() {
54 Matrix = NULL;
55 Indices = NULL;
56 Header = NULL;
57 MatrixCounter = 0;
58 RowCounter = NULL;
59 ColumnCounter = 0;
60};
61
62/** Destructor of MatrixContainer class.
63 */
64MatrixContainer::~MatrixContainer() {
65 if (Matrix != NULL)
66 for(int i=0;i<MatrixCounter;i++) {
67 if (RowCounter != NULL)
68 for(int j=0;j<=RowCounter[i];j++)
69 Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
70 Free((void **)&Matrix[i], "MatrixContainer::~MatrixContainer: **Matrix[]");
71 }
72 if ((RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
73 for(int j=0;j<=RowCounter[MatrixCounter];j++)
74 Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
75 if (MatrixCounter != 0)
76 Free((void **)&Matrix[MatrixCounter], "MatrixContainer::~MatrixContainer: **Matrix[MatrixCounter]");
77 Free((void **)&Matrix, "MatrixContainer::~MatrixContainer: ***Matrix");
78
79 if (Indices != NULL)
80 for(int i=0;i<=MatrixCounter;i++) {
81 Free((void **)&Indices[i], "MatrixContainer::~MatrixContainer: *Indices[]");
82 }
83 Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices");
84
85 Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
86 Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
87};
88
89/** Parsing a number of matrices.
90 * \param *name directory with files
91 * \param *prefix prefix of each matrix file
92 * \param *suffix suffix of each matrix file
93 * \param skiplines number of inital lines to skip
94 * \param skiplines number of inital columns to skip
95 * \return parsing successful
96 */
97bool MatrixContainer::ParseMatrix(char *name, char *prefix, char *suffix, int skiplines, int skipcolumns)
98{
99 char filename[1023];
100 ifstream input;
101 char *FragmentNumber = NULL;
102 stringstream line;
103 double tmp;
104
105 Header = (char *) Malloc(sizeof(char)*1023, "MatrixContainer::ParseMatrix: *EnergyHeader");
106
107 // count the number of matrices
108 MatrixCounter = -1; // we count one too much
109 line << name << FRAGMENTPREFIX << KEYSETFILE;
110 input.open(line.str().c_str(), ios::in);
111 if (input == NULL) {
112 cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
113 return false;
114 }
115 while (!input.eof()) {
116 input >> tmp;
117 MatrixCounter++;
118 }
119 input.close();
120 cout << "Determined " << MatrixCounter << " fragments." << endl;
121
122 cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
123 Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseMatrix: ***Matrix"); // one more each for the total molecule
124 RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseMatrix: *RowCounter");
125 for(int i=0;i<=MatrixCounter;i++) {
126 // energy
127 stringstream line;
128 FragmentNumber = FixedDigitNumber(MatrixCounter, i);
129 if (i != MatrixCounter)
130 line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix << suffix;
131 else
132 line << name << prefix << suffix;
133 Free((void **)&FragmentNumber, "MatrixContainer::ParseMatrix: *FragmentNumber");
134 input.open(line.str().c_str(), ios::in);
135 //cout << "Opening " << line.str() << " ... " << input << endl;
136 if (input == NULL) {
137 cerr << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
138 return false;
139 }
140 // scan header for number of columns
141 for (int m=0;m<=skiplines;m++)
142 input.getline(Header, 1023);
143 line.str(Header);
144 for(int k=0;k<skipcolumns;k++)
145 line >> Header;
146 //cout << line.str() << endl;
147 ColumnCounter=0;
148 while (!line.eof()) {
149 strcpy(filename,"@");
150 line >> filename;
151 if (filename[0] != '@')
152 ColumnCounter++;
153 }
154 //cout << line.str() << endl;
155 //cout << "ColumnCounter: " << ColumnCounter << endl;
156 // scan rest for number of rows/lines
157 RowCounter[i]=-1; // counts one line too much
158 while (!input.eof()) {
159 input.getline(filename, 1023);
160 //cout << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
161 RowCounter[i]++; // then line was not last MeanForce
162 if (strncmp(filename,"MeanForce", 9) == 0) {
163 break;
164 }
165 }
166 cout << "RowCounter[" << i << "]: " << RowCounter[i] << endl;
167 Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseMatrix: **Matrix[]");
168 input.clear();
169 input.seekg(ios::beg);
170 for (int m=0;m<=skiplines;m++)
171 input.getline(Header, 1023); // skip header
172 line.str(Header);
173 for(int k=0;k<skipcolumns;k++) // skip columns in header too
174 line >> filename;
175 strncpy(Header, line.str().c_str(), 1023);
176 for(int j=0;j<RowCounter[i];j++) {
177 Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseMatrix: *Matrix[][]");
178 input.getline(filename, 1023);
179 stringstream line(filename);
180 //cout << "Matrix at level " << j << ":";// << filename << endl;
181 for(int k=0;k<skipcolumns;k++)
182 line >> filename;
183 for(int k=0;(k<ColumnCounter) && (!line.eof());k++) {
184 line >> Matrix[i][j][k];
185 //cout << " " << setprecision(2) << Matrix[i][j][k];;
186 }
187 //cout << endl;
188 }
189 Matrix[i][ RowCounter[i] ] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseMatrix: *Matrix[RowCounter[i]][]");
190 for(int j=0;j<ColumnCounter;j++)
191 Matrix[i][ RowCounter[i] ][j] = 0.;
192 input.close();
193 }
194 return true;
195};
196
197/** Allocates and resets the memory for a number \a MCounter of matrices.
198 * \param *GivenHeader Header line
199 * \param MCounter number of matrices
200 * \param *RCounter number of rows for each matrix
201 * \param CCounter number of columns for all matrices
202 * \return Allocation successful
203 */
204bool MatrixContainer::AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int CCounter)
205{
206 Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseMatrix: *EnergyHeader");
207 strncpy(Header, GivenHeader, 1023);
208
209 MatrixCounter = MCounter;
210 ColumnCounter = CCounter;
211 Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseMatrix: ***Matrix"); // one more each for the total molecule
212 RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseMatrix: *RowCounter");
213 for(int i=0;i<=MatrixCounter;i++) {
214 RowCounter[i] = RCounter[i];
215 Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseMatrix: **Matrix[]");
216 for(int j=0;j<=RowCounter[i];j++) {
217 Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseMatrix: *Matrix[][]");
218 for(int k=0;k<ColumnCounter;k++)
219 Matrix[i][j][k] = 0.;
220 }
221 }
222 return true;
223};
224
225/** Resets all values in MatrixContainer::Matrix.
226 * \return true if successful
227 */
228bool MatrixContainer::ResetMatrix()
229{
230 for(int i=0;i<=MatrixCounter;i++)
231 for(int j=0;j<=RowCounter[i];j++)
232 for(int k=0;k<ColumnCounter;k++)
233 Matrix[i][j][k] = 0.;
234 return true;
235};
236
237/** Sets all values in the last of MatrixContainer::Matrix to \a value.
238 * \param value reset value
239 * \param skipcolumns skip initial columns
240 * \return true if successful
241 */
242bool MatrixContainer::SetLastMatrix(double value, int skipcolumns)
243{
244 for(int j=0;j<=RowCounter[MatrixCounter];j++)
245 for(int k=skipcolumns;k<ColumnCounter;k++)
246 Matrix[MatrixCounter][j][k] = value;
247 return true;
248};
249
250/** Sets all values in the last of MatrixContainer::Matrix to \a value.
251 * \param **values matrix with each value (must have at least same dimensions!)
252 * \param skipcolumns skip initial columns
253 * \return true if successful
254 */
255bool MatrixContainer::SetLastMatrix(double **values, int skipcolumns)
256{
257 for(int j=0;j<=RowCounter[MatrixCounter];j++)
258 for(int k=skipcolumns;k<ColumnCounter;k++)
259 Matrix[MatrixCounter][j][k] = values[j][k];
260 return true;
261};
262
263/** Sums the energy with each factor and put into last element of \a ***Energies.
264 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
265 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
266 * \param Order bond order
267 * \return true if summing was successful
268 */
269bool MatrixContainer::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
270{
271 // go through each order
272 for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
273 //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
274 // then go per order through each suborder and pick together all the terms that contain this fragment
275 for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
276 for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
277 if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
278 //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
279 // if the fragment's indices are all in the current fragment
280 for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
281 int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
282 //cout << "Current index is " << k << "/" << m << "." << endl;
283 if (m != -1) { // if it's not an added hydrogen
284 for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
285 //cout << "Comparing " << m << " with " << Matrix.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
286 if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
287 m = l;
288 break;
289 }
290 }
291 //cout << "Corresponding index in CurrentFragment is " << m << "." << endl;
292 if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
293 cerr << "Current force index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
294 return false;
295 }
296 if (Order == SubOrder) { // equal order is always copy from Energies
297 for(int l=0;l<ColumnCounter;l++) // then adds/subtract each column
298 Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
299 } else {
300 for(int l=0;l<ColumnCounter;l++)
301 Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
302 }
303 }
304 //if ((ColumnCounter>1) && (RowCounter >= 1))
305 //cout << "Fragments[ KeySet.OrderSet[" << i << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " << Fragments[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
306 }
307 } else {
308 //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet. KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
309 }
310 }
311 }
312 //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " << Fragments[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
313 }
314
315 return true;
316};
317
318/** Writes the summed total fragment terms \f$F_{ij}\f$ to file.
319 * \param *name inputdir
320 * \param *prefix prefix before \a EnergySuffix
321 * \return file was written
322 */
323bool MatrixContainer::WriteTotalFragments(const char *name, const char *prefix)
324{
325 ofstream output;
326 char *FragmentNumber = NULL;
327
328 cout << "Writing fragment files." << endl;
329 for(int i=0;i<MatrixCounter;i++) {
330 stringstream line;
331 FragmentNumber = FixedDigitNumber(MatrixCounter, i);
332 line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
333 Free((void **)&FragmentNumber, "*FragmentNumber");
334 output.open(line.str().c_str(), ios::out);
335 if (output == NULL) {
336 cerr << "Unable to open output energy file " << line.str() << "!" << endl;
337 return false;
338 }
339 output << Header << endl;
340 for(int j=0;j<RowCounter[i];j++) {
341 for(int k=0;k<ColumnCounter;k++)
342 output << scientific << Matrix[i][j][k] << "\t";
343 output << endl;
344 }
345 output.close();
346 }
347 return true;
348};
349
350/** Writes the summed total values in the last matrix to file.
351 * \param *name inputdir
352 * \param *prefix prefix
353 * \param *suffix suffix
354 * \return file was written
355 */
356bool MatrixContainer::WriteLastMatrix(const char *name, const char *prefix, const char *suffix)
357{
358 ofstream output;
359 stringstream line;
360
361 cout << "Writing matrix values of " << suffix << "." << endl;
362 line << name << prefix << suffix;
363 output.open(line.str().c_str(), ios::out);
364 if (output == NULL) {
365 cerr << "Unable to open output matrix file " << line.str() << "!" << endl;
366 return false;
367 }
368 output << Header << endl;
369 for(int j=0;j<RowCounter[MatrixCounter];j++) {
370 for(int k=0;k<ColumnCounter;k++)
371 output << scientific << Matrix[MatrixCounter][j][k] << "\t";
372 output << endl;
373 }
374 output.close();
375 return true;
376};
377
378// ======================================= CLASS EnergyMatrix =============================
379
380/** Create a trivial energy index mapping.
381 * This just maps 1 to 1, 2 to 2 and so on for all fragments.
382 * \return creation sucessful
383 */
384bool EnergyMatrix::ParseIndices()
385{
386 cout << "Parsing energy indices." << endl;
387 Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "EnergyMatrix::ParseIndices: **Indices");
388 for(int i=0;i<=MatrixCounter;i++) {
389 Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "EnergyMatrix::ParseIndices: *Indices[]");
390 for(int j=0;j<RowCounter[i];j++)
391 Indices[i][j] = j;
392 }
393 return true;
394};
395
396/** Sums the energy with each factor and put into last element of \a EnergyMatrix::Matrix.
397 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
398 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
399 * \param Order bond order
400 * \return true if summing was successful
401 */
402bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, int Order)
403{
404 // sum energy
405 for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++)
406 for(int j=0;j<RowCounter[ KeySet.OrderSet[Order][i] ];j++)
407 for(int k=0;k<ColumnCounter;k++)
408 Matrix[MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
409 return true;
410};
411
412// ======================================= CLASS ForceMatrix =============================
413
414/** Parsing force Indices of each fragment
415 * \param *name directory with \a ForcesFile
416 * \return parsing successful
417 */
418bool ForceMatrix::ParseIndices(char *name)
419{
420 cout << "Parsing force indices." << endl;
421 Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "ForceMatrix::ParseIndices: **Indices");
422 for (int i=0;i<MatrixCounter;i++) {
423 Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
424 for(int j=0;j<RowCounter[i];j++)
425 Indices[i][j] = j;
426 }
427 Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
428 for(int j=0;j<RowCounter[MatrixCounter];j++) {
429 Indices[MatrixCounter][j] = j;
430 }
431 return true;
432};
433
434
435/** Sums the forces and puts into last element of \a ForceMatrix::Matrix.
436 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
437 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
438 * \param Order bond order
439 * \return true if summing was successful
440 */
441bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order)
442{
443 // sum forces
444 for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++)
445 for(int l=0;l<RowCounter[ KeySet.OrderSet[Order][i] ];l++) {
446 int j = Indices[ KeySet.OrderSet[Order][i] ][l];
447 if (j > RowCounter[MatrixCounter]) {
448 cerr << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
449 return false;
450 }
451 if (j != -1)
452 for(int k=2;k<ColumnCounter;k++)
453 Matrix[MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[Order][i] ][l][k];
454 }
455 return true;
456};
457
458
459// ======================================= CLASS KeySetsContainer =============================
460
461/** Constructor of KeySetsContainer class.
462 */
463KeySetsContainer::KeySetsContainer() {
464 KeySets = NULL;
465 AtomCounter = NULL;
466 FragmentCounter = 0;
467 Order = 0;
468 FragmentsPerOrder = 0;
469 OrderSet = NULL;
470};
471
472/** Destructor of KeySetsContainer class.
473 */
474KeySetsContainer::~KeySetsContainer() {
475 for(int i=0;i<FragmentCounter;i++)
476 Free((void **)&KeySets[i], "KeySetsContainer::~KeySetsContainer: *KeySets[]");
477 for(int i=0;i<Order;i++)
478 Free((void **)&OrderSet[i], "KeySetsContainer::~KeySetsContainer: *OrderSet[]");
479 Free((void **)&KeySets, "KeySetsContainer::~KeySetsContainer: **KeySets");
480 Free((void **)&OrderSet, "KeySetsContainer::~KeySetsContainer: **OrderSet");
481 Free((void **)&AtomCounter, "KeySetsContainer::~KeySetsContainer: *AtomCounter");
482 Free((void **)&FragmentsPerOrder, "KeySetsContainer::~KeySetsContainer: *FragmentsPerOrder");
483};
484
485/** Parsing KeySets into array.
486 * \param *name directory with keyset file
487 * \param *ACounter number of atoms per fragment
488 * \param FCounter number of fragments
489 * \return parsing succesful
490 */
491bool KeySetsContainer::ParseKeySets(const char *name, const int *ACounter, const int FCounter) {
492 ifstream input;
493 char *FragmentNumber = NULL;
494 stringstream line;
495 char filename[1023];
496
497 FragmentCounter = FCounter;
498 cout << "Parsing key sets." << endl;
499 KeySets = (int **) Malloc(sizeof(int *)*FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
500 line << name << FRAGMENTPREFIX << KEYSETFILE;
501 input.open(line.str().c_str(), ios::in);
502 if (input == NULL) {
503 cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
504 return false;
505 }
506
507 AtomCounter = (int *) Malloc(sizeof(int)*FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
508 for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
509 stringstream line;
510 AtomCounter[i] = ACounter[i];
511 // parse the values
512 KeySets[i] = (int *) Malloc(sizeof(int)*AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
513 FragmentNumber = FixedDigitNumber(FragmentCounter, i);
514 cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
515 Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
516 input.getline(filename, 1023);
517 line.str(filename);
518 for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
519 line >> KeySets[i][j];
520 cout << " " << KeySets[i][j];
521 }
522 cout << endl;
523 }
524 input.close();
525 return true;
526};
527
528/** Parse many body terms, associating each fragment to a certain bond order.
529 * \return parsing succesful
530 */
531bool KeySetsContainer::ParseManyBodyTerms()
532{
533 int Counter;
534
535 cout << "Creating Fragment terms." << endl;
536 // scan through all to determine order
537 Order=0;
538 for(int i=0;i<FragmentCounter;i++) {
539 Counter=0;
540 for(int j=0;j<AtomCounter[i];j++)
541 if (KeySets[i][j] != -1)
542 Counter++;
543 if (Counter > Order)
544 Order = Counter;
545 }
546 cout << "Found Order is " << Order << "." << endl;
547
548 // scan through all to determine fragments per order
549 FragmentsPerOrder = (int *) Malloc(sizeof(int)*Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
550 for(int i=0;i<Order;i++)
551 FragmentsPerOrder[i] = 0;
552 for(int i=0;i<FragmentCounter;i++) {
553 Counter=0;
554 for(int j=0;j<AtomCounter[i];j++)
555 if (KeySets[i][j] != -1)
556 Counter++;
557 FragmentsPerOrder[Counter-1]++;
558 }
559 for(int i=0;i<Order;i++)
560 cout << "Found No. of Fragments of Order " << i+1 << " is " << FragmentsPerOrder[i] << "." << endl;
561
562 // scan through all to gather indices to each order set
563 OrderSet = (int **) Malloc(sizeof(int *)*Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
564 for(int i=0;i<Order;i++)
565 OrderSet[i] = (int *) Malloc(sizeof(int)*FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
566 for(int i=0;i<Order;i++)
567 FragmentsPerOrder[i] = 0;
568 for(int i=0;i<FragmentCounter;i++) {
569 Counter=0;
570 for(int j=0;j<AtomCounter[i];j++)
571 if (KeySets[i][j] != -1)
572 Counter++;
573 OrderSet[Counter-1][FragmentsPerOrder[Counter-1]] = i;
574 FragmentsPerOrder[Counter-1]++;
575 }
576 cout << "Printing OrderSet." << endl;
577 for(int i=0;i<Order;i++) {
578 for (int j=0;j<FragmentsPerOrder[i];j++) {
579 cout << " " << OrderSet[i][j];
580 }
581 cout << endl;
582 }
583 cout << endl;
584
585
586 return true;
587};
588
589/** Compares each entry in \a *SmallerSet if it is containted in \a *GreaterSet.
590 * \param GreaterSet index to greater set
591 * \param SmallerSet index to smaller set
592 * \return true if all keys of SmallerSet contained in GreaterSet
593 */
594bool KeySetsContainer::Contains(const int GreaterSet, const int SmallerSet)
595{
596 bool result = true;
597 bool intermediate;
598 if ((GreaterSet < 0) || (SmallerSet < 0) || (GreaterSet > FragmentCounter) || (SmallerSet > FragmentCounter)) // index out of bounds
599 return false;
600 for(int i=0;i<AtomCounter[SmallerSet];i++) {
601 intermediate = false;
602 for (int j=0;j<AtomCounter[GreaterSet];j++)
603 intermediate = (intermediate || ((KeySets[SmallerSet][i] == KeySets[GreaterSet][j]) || (KeySets[SmallerSet][i] == -1)));
604 result = result && intermediate;
605 }
606
607 return result;
608};
609
610
611// ======================================= END =============================================
Note: See TracBrowser for help on using the repository browser.