source: src/parser.cpp@ 6f6a8e

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

some comments for where "F" and where "E" terms are created

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