source: src/parser.cpp@ 2459b1

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 2459b1 was 2459b1, checked in by Frederik Heber <heber@…>, 17 years ago

ForcesFile and TEFactors are _needed_, reincorporated. UniqueFragments structure is now in BOSSANOVA

ForcesFile is again written, as we need it for the hydrogen not coming saturation (forces!)
TEFactors are back, as despite my notion they are needed in the evaluation
UniqueFragments structure is shifted from PowerSetGenerator to FragmentBOSSANOVA. Actually only for the labels - however, the if was changed to test the real indices (which are also always the same, which is better for adaptive runs) - but might be more useful there still.
analyzer and joiner again parse indices.

  • Property mode set to 100644
File size: 23.2 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 << KEYSETFILE;
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.getline(filename, 1023);
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.