source: src/moleculelist.cpp@ fae462

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 fae462 was 910a5d, checked in by Frederik Heber <heber@…>, 10 years ago

atom::father is now accessible through getter and setter only.

  • this allows to manage signing on and off to and from father to be notified when it's destroyed.
  • Property mode set to 100755
File size: 21.9 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
5 *
6 *
7 * This file is part of MoleCuilder.
8 *
9 * MoleCuilder is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * MoleCuilder is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
21 */
22
23/** \file MoleculeListClass.cpp
24 *
25 * Function implementations for the class MoleculeListClass.
26 *
27 */
28
29// include config.h
30#ifdef HAVE_CONFIG_H
31#include <config.h>
32#endif
33
34#include "CodePatterns/MemDebug.hpp"
35
36#include <iostream>
37
38//#include <gsl/gsl_inline.h>
39#include <gsl/gsl_heapsort.h>
40
41#include "MoleculeListClass.hpp"
42
43#include "CodePatterns/Log.hpp"
44
45#include "Atom/atom.hpp"
46#include "Bond/bond.hpp"
47#include "Box.hpp"
48#include "config.hpp"
49#include "Element/element.hpp"
50#include "Element/periodentafel.hpp"
51#include "Fragmentation/Graph.hpp"
52#include "Fragmentation/KeySet.hpp"
53#include "Fragmentation/SortIndex.hpp"
54#include "Graph/BondGraph.hpp"
55#include "Helpers/helpers.hpp"
56#include "molecule.hpp"
57#include "LinearAlgebra/RealSpaceMatrix.hpp"
58#include "Parser/FormatParserStorage.hpp"
59#include "World.hpp"
60
61
62/** Constructor for MoleculeListClass.
63 */
64MoleculeListClass::MoleculeListClass(World *_world) :
65 Observable("MoleculeListClass"),
66 MaxIndex(1),
67 world(_world)
68{};
69
70/** Destructor for MoleculeListClass.
71 */
72MoleculeListClass::~MoleculeListClass()
73{
74 LOG(4, "Clearing ListOfMolecules.");
75 for(MoleculeList::iterator MolRunner = ListOfMolecules.begin(); MolRunner != ListOfMolecules.end(); ++MolRunner)
76 (*MolRunner)->signOff(this);
77 ListOfMolecules.clear(); // empty list
78};
79
80/** Insert a new molecule into the list and set its number.
81 * \param *mol molecule to add to list.
82 */
83void MoleculeListClass::insert(molecule *mol)
84{
85 OBSERVE;
86 mol->IndexNr = MaxIndex++;
87 ListOfMolecules.push_back(mol);
88 mol->signOn(this);
89};
90
91/** Erases a molecule from the list.
92 * \param *mol molecule to add to list.
93 */
94void MoleculeListClass::erase(molecule *mol)
95{
96 OBSERVE;
97 mol->signOff(this);
98 ListOfMolecules.remove(mol);
99};
100
101/** Comparison function for two values.
102 * \param *a
103 * \param *b
104 * \return <0, \a *a less than \a *b, ==0 if equal, >0 \a *a greater than \a *b
105 */
106int CompareDoubles (const void * a, const void * b)
107{
108 if (*(double *)a > *(double *)b)
109 return -1;
110 else if (*(double *)a < *(double *)b)
111 return 1;
112 else
113 return 0;
114};
115
116
117/** Compare whether two molecules are equal.
118 * \param *a molecule one
119 * \param *n molecule two
120 * \return lexical value (-1, 0, +1)
121 */
122int MolCompare(const void *a, const void *b)
123{
124 int *aList = NULL, *bList = NULL;
125 int Count, Counter, aCounter, bCounter;
126 int flag;
127
128 // sort each atom list and put the numbers into a list, then go through
129 //LOG(0, "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << ".");
130 // Yes those types are awkward... but check it for yourself it checks out this way
131 molecule *const *mol1_ptr= static_cast<molecule *const *>(a);
132 molecule *mol1 = *mol1_ptr;
133 molecule *const *mol2_ptr= static_cast<molecule *const *>(b);
134 molecule *mol2 = *mol2_ptr;
135 if (mol1->getAtomCount() < mol2->getAtomCount()) {
136 return -1;
137 } else {
138 if (mol1->getAtomCount() > mol2->getAtomCount())
139 return +1;
140 else {
141 Count = mol1->getAtomCount();
142 aList = new int[Count];
143 bList = new int[Count];
144
145 // fill the lists
146 Counter = 0;
147 aCounter = 0;
148 bCounter = 0;
149 molecule::const_iterator aiter = mol1->begin();
150 molecule::const_iterator biter = mol2->begin();
151 for (;(aiter != mol1->end()) && (biter != mol2->end());
152 ++aiter, ++biter) {
153 if ((*aiter)->GetTrueFather() == NULL)
154 aList[Counter] = Count + (aCounter++);
155 else
156 aList[Counter] = (*aiter)->GetTrueFather()->getNr();
157 if ((*biter)->GetTrueFather() == NULL)
158 bList[Counter] = Count + (bCounter++);
159 else
160 bList[Counter] = (*biter)->GetTrueFather()->getNr();
161 Counter++;
162 }
163 // check if AtomCount was for real
164 flag = 0;
165 if ((aiter == mol1->end()) && (biter != mol2->end())) {
166 flag = -1;
167 } else {
168 if ((aiter != mol1->end()) && (biter == mol2->end()))
169 flag = 1;
170 }
171 if (flag == 0) {
172 // sort the lists
173 gsl_heapsort(aList, Count, sizeof(int), CompareDoubles);
174 gsl_heapsort(bList, Count, sizeof(int), CompareDoubles);
175 // compare the lists
176
177 flag = 0;
178 for (int i = 0; i < Count; i++) {
179 if (aList[i] < bList[i]) {
180 flag = -1;
181 } else {
182 if (aList[i] > bList[i])
183 flag = 1;
184 }
185 if (flag != 0)
186 break;
187 }
188 }
189 delete[] (aList);
190 delete[] (bList);
191 return flag;
192 }
193 }
194 return -1;
195};
196
197/** Output of a list of all molecules.
198 * \param *out output stream
199 */
200void MoleculeListClass::Enumerate(std::ostream *out)
201{
202 periodentafel *periode = World::getInstance().getPeriode();
203 std::map<atomicNumber_t,unsigned int> counts;
204 double size=0;
205 Vector Origin;
206
207 // header
208 (*out) << "Index\tName\t\tAtoms\tFormula\tCenter\tSize" << endl;
209 (*out) << "-----------------------------------------------" << endl;
210 if (ListOfMolecules.size() == 0)
211 (*out) << "\tNone" << endl;
212 else {
213 Origin.Zero();
214 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
215 // count atoms per element and determine size of bounding sphere
216 size=0.;
217 for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
218 counts[(*iter)->getType()->getAtomicNumber()]++;
219 if ((*iter)->DistanceSquared(Origin) > size)
220 size = (*iter)->DistanceSquared(Origin);
221 }
222 // output Index, Name, number of atoms, chemical formula
223 (*out) << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->getAtomCount() << "\t";
224
225 std::map<atomicNumber_t,unsigned int>::reverse_iterator iter;
226 for(iter=counts.rbegin(); iter!=counts.rend();++iter){
227 atomicNumber_t Z =(*iter).first;
228 (*out) << periode->FindElement(Z)->getSymbol() << (*iter).second;
229 }
230 // Center and size
231 Vector Center = (*ListRunner)->DetermineCenterOfAll();
232 (*out) << "\t" << Center << "\t" << sqrt(size) << endl;
233 }
234 }
235};
236
237/** Returns the molecule with the given index \a index.
238 * \param index index of the desired molecule
239 * \return pointer to molecule structure, NULL if not found
240 */
241molecule * MoleculeListClass::ReturnIndex(int index)
242{
243 for(MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
244 if ((*ListRunner)->IndexNr == index)
245 return (*ListRunner);
246 return NULL;
247};
248
249
250/** Simple output of the pointers in ListOfMolecules.
251 * \param *out output stream
252 */
253void MoleculeListClass::Output(std::ostream *out)
254{
255 if (DoLog(1)) {
256 std::stringstream output;
257 output << "MoleculeList: ";
258 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
259 output << *ListRunner << "\t";
260 LOG(1, output.str());
261 }
262};
263
264/** Returns a string with \a i prefixed with 0s to match order of total number of molecules in digits.
265 * \param FragmentNumber total number of fragments to determine necessary number of digits
266 * \param digits number to create with 0 prefixed
267 * \return allocated(!) char array with number in digits, ten base.
268 */
269inline char *FixedDigitNumber(const int FragmentNumber, const int digits)
270{
271 char *returnstring;
272 int number = FragmentNumber;
273 int order = 0;
274 while (number != 0) { // determine number of digits needed
275 number = (int)floor(((double)number / 10.));
276 order++;
277 //LOG(0, "Number is " << number << ", order is " << order << ".");
278 }
279 // allocate string
280 returnstring = new char[order + 2];
281 // terminate and fill string array from end backward
282 returnstring[order] = '\0';
283 number = digits;
284 for (int i=order;i--;) {
285 returnstring[i] = '0' + (char)(number % 10);
286 number = (int)floor(((double)number / 10.));
287 }
288 //LOG(0, returnstring);
289 return returnstring;
290};
291
292/** Calculates necessary hydrogen correction due to unwanted interaction between saturated ones.
293 * If for a pair of two hydrogen atoms a and b, at least is a saturated one, and a and b are not
294 * bonded to the same atom, then we add for this pair a correction term constructed from a Morse
295 * potential function fit to QM calculations with respecting to the interatomic hydrogen distance.
296 * \param &path path to file
297 */
298bool MoleculeListClass::AddHydrogenCorrection(std::string &path)
299{
300 bond::ptr Binder;
301 double ***FitConstant = NULL, **correction = NULL;
302 int a, b;
303 ofstream output;
304 ifstream input;
305 string line;
306 stringstream zeile;
307 double distance;
308 char ParsedLine[1023];
309 double tmp;
310 char *FragmentNumber = NULL;
311
312 LOG(1, "Saving hydrogen saturation correction ... ");
313 // 0. parse in fit constant files that should have the same dimension as the final energy files
314 // 0a. find dimension of matrices with constants
315 line = path;
316 line += "1";
317 line += FITCONSTANTSUFFIX;
318 input.open(line.c_str());
319 if (input.fail()) {
320 LOG(1, endl << "Unable to open " << line << ", is the directory correct?");
321 return false;
322 }
323 a = 0;
324 b = -1; // we overcount by one
325 while (!input.eof()) {
326 input.getline(ParsedLine, 1023);
327 zeile.str(ParsedLine);
328 int i = 0;
329 while (!zeile.eof()) {
330 zeile >> distance;
331 i++;
332 }
333 if (i > a)
334 a = i;
335 b++;
336 }
337 LOG(0, "I recognized " << a << " columns and " << b << " rows, ");
338 input.close();
339
340 // 0b. allocate memory for constants
341 FitConstant = new double**[3];
342 for (int k = 0; k < 3; k++) {
343 FitConstant[k] = new double*[a];
344 for (int i = a; i--;) {
345 FitConstant[k][i] = new double[b];
346 for (int j = b; j--;) {
347 FitConstant[k][i][j] = 0.;
348 }
349 }
350 }
351 // 0c. parse in constants
352 for (int i = 0; i < 3; i++) {
353 line = path;
354 line.append("/");
355 line += FRAGMENTPREFIX;
356 sprintf(ParsedLine, "%d", i + 1);
357 line += ParsedLine;
358 line += FITCONSTANTSUFFIX;
359 input.open(line.c_str());
360 if (input.fail()) {
361 ELOG(0, endl << "Unable to open " << line << ", is the directory correct?");
362 performCriticalExit();
363 return false;
364 }
365 int k = 0, l;
366 while ((!input.eof()) && (k < b)) {
367 input.getline(ParsedLine, 1023);
368 //LOG(1, "INFO: Current Line: " << ParsedLine);
369 zeile.str(ParsedLine);
370 zeile.clear();
371 l = 0;
372 //std::stringstream output;
373 while ((!zeile.eof()) && (l < a)) {
374 zeile >> FitConstant[i][l][k];
375 //output << FitConstant[i][l][k] << "\t";
376 l++;
377 }
378 //LOG(1, "INFO: fit constant are " << output.str());
379 k++;
380 }
381 input.close();
382 }
383 if (DoLog(1)) {
384 for (int k = 0; k < 3; k++) {
385 std::stringstream output;
386 output << "Constants " << k << ": ";
387 for (int j = 0; j < b; j++) {
388 for (int i = 0; i < a; i++) {
389 output << FitConstant[k][i][j] << "\t";
390 }
391 output << std::endl;
392 }
393 LOG(0, output.str());
394 }
395 }
396
397 // 0d. allocate final correction matrix
398 correction = new double*[a];
399 for (int i = a; i--;)
400 correction[i] = new double[b];
401
402 // 1a. go through every molecule in the list
403 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
404 // 1b. zero final correction matrix
405 for (int k = a; k--;)
406 for (int j = b; j--;)
407 correction[k][j] = 0.;
408 // 2. take every hydrogen that is a saturated one
409 for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
410 //LOG(1, "(*iter): " << *(*iter) << " with first bond " << *((*iter)->getListOfBonds().begin()) << ".");
411 if (((*iter)->getType()->getAtomicNumber() == 1) && (((*iter)->getFather() == NULL)
412 || ((*iter)->getFather()->getType()->getAtomicNumber() != 1))) { // if it's a hydrogen
413 for (molecule::const_iterator runner = (*ListRunner)->begin(); runner != (*ListRunner)->end(); ++runner) {
414 //LOG(2, "Runner: " << *(*runner) << " with first bond " << *((*iter)->getListOfBonds().begin()) << ".");
415 // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
416 const BondList &bondlist = (*runner)->getListOfBonds();
417 Binder = *(bondlist.begin());
418 if (((*runner)->getType()->getAtomicNumber() == 1) && ((*runner)->getNr() > (*iter)->getNr()) && (Binder->GetOtherAtom((*runner)) != Binder->GetOtherAtom((*iter)))) { // (hydrogens have only one bonding partner!)
419 // 4. evaluate the morse potential for each matrix component and add up
420 distance = (*runner)->distance(*(*iter));
421 //std::stringstream output;
422 //output << "Fragment " << (*ListRunner)->name << ": " << *(*runner) << "<= " << distance << "=>" << *(*iter) << ":";
423 for (int k = 0; k < a; k++) {
424 for (int j = 0; j < b; j++) {
425 switch (k) {
426 case 1:
427 case 7:
428 case 11:
429 tmp = pow(FitConstant[0][k][j] * (1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]))), 2);
430 break;
431 default:
432 tmp = FitConstant[0][k][j] * pow(distance, FitConstant[1][k][j]) + FitConstant[2][k][j];
433 break;
434 };
435 correction[k][j] -= tmp; // ground state is actually lower (disturbed by additional interaction)
436 //output << tmp << "\t";
437 }
438 //output << endl;
439 }
440 //LOG(0, output.str());
441 }
442 }
443 }
444 }
445 // 5. write final matrix to file
446 line = path;
447 line.append("/");
448 line += FRAGMENTPREFIX;
449 FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), (*ListRunner)->IndexNr);
450 line += FragmentNumber;
451 delete[] (FragmentNumber);
452 line += HCORRECTIONSUFFIX;
453 output.open(line.c_str());
454 output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
455 for (int j = 0; j < b; j++) {
456 for (int i = 0; i < a; i++)
457 output << correction[i][j] << "\t";
458 output << endl;
459 }
460 output.close();
461 }
462 for (int i = a; i--;)
463 delete[](correction[i]);
464 delete[](correction);
465
466 line = path;
467 line.append("/");
468 line += HCORRECTIONSUFFIX;
469 output.open(line.c_str());
470 output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
471 for (int j = 0; j < b; j++) {
472 for (int i = 0; i < a; i++)
473 output << 0 << "\t";
474 output << endl;
475 }
476 output.close();
477 // 6. free memory of parsed matrices
478 for (int k = 0; k < 3; k++) {
479 for (int i = a; i--;) {
480 delete[](FitConstant[k][i]);
481 }
482 delete[](FitConstant[k]);
483 }
484 delete[](FitConstant);
485 LOG(0, "done.");
486 return true;
487};
488
489/** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
490 * \param &path path to file
491 * \param SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
492 * \return true - file written successfully, false - writing failed
493 */
494bool MoleculeListClass::StoreForcesFile(std::string &path, const SortIndex_t &SortIndex)
495{
496 bool status = true;
497 string filename(path);
498 filename += FORCESFILE;
499 ofstream ForcesFile(filename.c_str());
500 periodentafel *periode=World::getInstance().getPeriode();
501
502 // open file for the force factors
503 LOG(1, "Saving force factors ... ");
504 if (!ForcesFile.fail()) {
505 //output << prefix << "Forces" << endl;
506 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
507 periodentafel::const_iterator elemIter;
508 for(elemIter=periode->begin();elemIter!=periode->end();++elemIter){
509 if ((*ListRunner)->hasElement((*elemIter).first)) { // if this element got atoms
510 for(molecule::iterator atomIter = (*ListRunner)->begin(); atomIter !=(*ListRunner)->end();++atomIter){
511 if ((*atomIter)->getType()->getAtomicNumber() == (*elemIter).first) {
512 if (((*atomIter)->GetTrueFather() != NULL) && ((*atomIter)->GetTrueFather() != (*atomIter))) {// if there is a rea
513 const atomId_t fatherid = (*atomIter)->GetTrueFather()->getId();
514 ForcesFile << SortIndex.find(fatherid) << "\t";
515 } else
516 // otherwise a -1 to indicate an added saturation hydrogen
517 ForcesFile << "-1\t";
518 }
519 }
520 }
521 }
522 ForcesFile << endl;
523 }
524 ForcesFile.close();
525 LOG(1, "done.");
526 } else {
527 status = false;
528 LOG(1, "failed to open file " << filename << ".");
529 }
530 ForcesFile.close();
531
532 return status;
533};
534
535/** Writes a config file for each molecule in the given \a **FragmentList.
536 * \param *out output stream for debugging
537 * \param &prefix path and prefix to the fragment config files
538 * \param SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
539 * \param type desired type to store
540 * \return true - success (each file was written), false - something went wrong.
541 */
542bool MoleculeListClass::OutputConfigForListOfFragments(std::string &prefix, ParserTypes type)
543{
544 ofstream outputFragment;
545 std::string FragmentName;
546 bool result = true;
547 bool intermediateResult = true;
548 Vector BoxDimension;
549 char *FragmentNumber = NULL;
550 int FragmentCounter = 0;
551 ofstream output;
552 RealSpaceMatrix cell_size = World::getInstance().getDomain().getM();
553 RealSpaceMatrix cell_size_backup = cell_size;
554 int count=0;
555
556 // store the fragments as config and as xyz
557 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
558 // correct periodic
559 if ((*ListRunner)->ScanForPeriodicCorrection()) {
560 count++;
561 }
562
563 {
564 // list atoms in fragment for debugging
565 std::stringstream output;
566 output << "Contained atoms: ";
567 for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
568 output << (*iter)->getName() << " ";
569 }
570 LOG(2, output.str());
571 }
572
573 {
574// // center on edge
575// (*ListRunner)->CenterEdge(&BoxDimension);
576// for (int k = 0; k < NDIM; k++) // if one edge is to small, set at least to 1 angstroem
577// if (BoxDimension[k] < 1.)
578// BoxDimension[k] += 1.;
579// {
580// RealSpaceMatrix domain;
581// for(int i =0; i<NDIM;++i)
582// domain.at(i,i) = BoxDimension[i];
583// World::getInstance().setDomain(domain);
584// }
585// for (int k = 0; k < NDIM; k++) {
586// BoxDimension[k] = 2.5 * (World::getInstance().getConfig()->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
587// cell_size.at(k,k) = BoxDimension[k] * 2.;
588// }
589// World::getInstance().setDomain(cell_size);
590// (*ListRunner)->Translate(&BoxDimension);
591
592 // output file
593 std::vector<atom *> atoms;
594 // TODO: Convert iterator to const_iterator when FormatParserStorage::save() has vector<const atom *>
595 // We need iterator here because FormatParserStorage::save() need vector<atom *> not const refs.
596 for (molecule::iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
597 atoms.push_back(*iter);
598 }
599 FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), FragmentCounter++);
600 FragmentName = prefix + FragmentNumber + "." + FormatParserStorage::getInstance().getSuffixFromType(type);
601 outputFragment.open(FragmentName.c_str(), ios::out);
602 std::stringstream output;
603 output << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as XYZ ... ";
604 if ((intermediateResult = FormatParserStorage::getInstance().save(
605 outputFragment,
606 FormatParserStorage::getInstance().getSuffixFromType(type),
607 atoms)))
608 output << " done.";
609 else
610 output << " failed.";
611 LOG(3, output.str());
612 delete[](FragmentNumber);
613
614 result = result && intermediateResult;
615 outputFragment.close();
616 outputFragment.clear();
617 }
618 }
619 LOG(0, "STATUS: done.");
620
621 // printing final number
622 LOG(2, "INFO: Final number of fragments: " << FragmentCounter << ".");
623
624 // printing final number
625 LOG(2, "INFO: For " << count << " fragments periodic correction would have been necessary.");
626
627 // restore cell_size
628 World::getInstance().setDomain(cell_size_backup);
629
630 return result;
631};
632
633/** Counts the number of molecules with the molecule::ActiveFlag set.
634 * \return number of molecules with ActiveFlag set to true.
635 */
636int MoleculeListClass::NumberOfActiveMolecules()
637{
638 int count = 0;
639 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
640 count += ((*ListRunner)->ActiveFlag ? 1 : 0);
641 return count;
642};
643
644/** Count all atoms in each molecule.
645 * \return number of atoms in the MoleculeListClass.
646 * TODO: the inner loop should be done by some (double molecule::CountAtom()) function
647 */
648int MoleculeListClass::CountAllAtoms() const
649{
650 int AtomNo = 0;
651 for (MoleculeList::const_iterator MolWalker = ListOfMolecules.begin(); MolWalker != ListOfMolecules.end(); MolWalker++) {
652 AtomNo += (*MolWalker)->size();
653 }
654 return AtomNo;
655}
Note: See TracBrowser for help on using the repository browser.