source: src/moleculelist.cpp@ 9f9b5a

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 9f9b5a was cca9ef, checked in by Frederik Heber <heber@…>, 14 years ago

Renamed Matrix to RealSpaceMatrix.

  • class Matrix only represents 3x3 matrices, whereas we are now going to extend the class MatrixContent to contain arbritrary dimensions.
  • renamed class and file
  • Property mode set to 100755
File size: 36.1 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010 University of Bonn. All rights reserved.
5 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
6 */
7
8/** \file MoleculeListClass.cpp
9 *
10 * Function implementations for the class MoleculeListClass.
11 *
12 */
13
14// include config.h
15#ifdef HAVE_CONFIG_H
16#include <config.h>
17#endif
18
19#include "Helpers/MemDebug.hpp"
20
21#include <cstring>
22
23#include <gsl/gsl_inline.h>
24#include <gsl/gsl_heapsort.h>
25
26#include "World.hpp"
27#include "atom.hpp"
28#include "bond.hpp"
29#include "bondgraph.hpp"
30#include "boundary.hpp"
31#include "config.hpp"
32#include "element.hpp"
33#include "Helpers/helpers.hpp"
34#include "linkedcell.hpp"
35#include "lists.hpp"
36#include "Helpers/Verbose.hpp"
37#include "Helpers/Log.hpp"
38#include "molecule.hpp"
39#include "periodentafel.hpp"
40#include "tesselation.hpp"
41#include "Helpers/Assert.hpp"
42#include "LinearAlgebra/RealSpaceMatrix.hpp"
43#include "Box.hpp"
44#include "stackclass.hpp"
45
46#include "Helpers/Assert.hpp"
47
48/*********************************** Functions for class MoleculeListClass *************************/
49
50/** Constructor for MoleculeListClass.
51 */
52MoleculeListClass::MoleculeListClass(World *_world) :
53 Observable("MoleculeListClass"),
54 MaxIndex(1),
55 world(_world)
56{};
57
58/** Destructor for MoleculeListClass.
59 */
60MoleculeListClass::~MoleculeListClass()
61{
62 DoLog(4) && (Log() << Verbose(4) << "Clearing ListOfMolecules." << endl);
63 for(MoleculeList::iterator MolRunner = ListOfMolecules.begin(); MolRunner != ListOfMolecules.end(); ++MolRunner)
64 (*MolRunner)->signOff(this);
65 ListOfMolecules.clear(); // empty list
66};
67
68/** Insert a new molecule into the list and set its number.
69 * \param *mol molecule to add to list.
70 */
71void MoleculeListClass::insert(molecule *mol)
72{
73 OBSERVE;
74 mol->IndexNr = MaxIndex++;
75 ListOfMolecules.push_back(mol);
76 mol->signOn(this);
77};
78
79/** Erases a molecule from the list.
80 * \param *mol molecule to add to list.
81 */
82void MoleculeListClass::erase(molecule *mol)
83{
84 OBSERVE;
85 mol->signOff(this);
86 ListOfMolecules.remove(mol);
87};
88
89/** Compare whether two molecules are equal.
90 * \param *a molecule one
91 * \param *n molecule two
92 * \return lexical value (-1, 0, +1)
93 */
94int MolCompare(const void *a, const void *b)
95{
96 int *aList = NULL, *bList = NULL;
97 int Count, Counter, aCounter, bCounter;
98 int flag;
99
100 // sort each atom list and put the numbers into a list, then go through
101 //Log() << Verbose(0) << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
102 // Yes those types are awkward... but check it for yourself it checks out this way
103 molecule *const *mol1_ptr= static_cast<molecule *const *>(a);
104 molecule *mol1 = *mol1_ptr;
105 molecule *const *mol2_ptr= static_cast<molecule *const *>(b);
106 molecule *mol2 = *mol2_ptr;
107 if (mol1->getAtomCount() < mol2->getAtomCount()) {
108 return -1;
109 } else {
110 if (mol1->getAtomCount() > mol2->getAtomCount())
111 return +1;
112 else {
113 Count = mol1->getAtomCount();
114 aList = new int[Count];
115 bList = new int[Count];
116
117 // fill the lists
118 Counter = 0;
119 aCounter = 0;
120 bCounter = 0;
121 molecule::const_iterator aiter = mol1->begin();
122 molecule::const_iterator biter = mol2->begin();
123 for (;(aiter != mol1->end()) && (biter != mol2->end());
124 ++aiter, ++biter) {
125 if ((*aiter)->GetTrueFather() == NULL)
126 aList[Counter] = Count + (aCounter++);
127 else
128 aList[Counter] = (*aiter)->GetTrueFather()->nr;
129 if ((*biter)->GetTrueFather() == NULL)
130 bList[Counter] = Count + (bCounter++);
131 else
132 bList[Counter] = (*biter)->GetTrueFather()->nr;
133 Counter++;
134 }
135 // check if AtomCount was for real
136 flag = 0;
137 if ((aiter == mol1->end()) && (biter != mol2->end())) {
138 flag = -1;
139 } else {
140 if ((aiter != mol1->end()) && (biter == mol2->end()))
141 flag = 1;
142 }
143 if (flag == 0) {
144 // sort the lists
145 gsl_heapsort(aList, Count, sizeof(int), CompareDoubles);
146 gsl_heapsort(bList, Count, sizeof(int), CompareDoubles);
147 // compare the lists
148
149 flag = 0;
150 for (int i = 0; i < Count; i++) {
151 if (aList[i] < bList[i]) {
152 flag = -1;
153 } else {
154 if (aList[i] > bList[i])
155 flag = 1;
156 }
157 if (flag != 0)
158 break;
159 }
160 }
161 delete[] (aList);
162 delete[] (bList);
163 return flag;
164 }
165 }
166 return -1;
167};
168
169/** Output of a list of all molecules.
170 * \param *out output stream
171 */
172void MoleculeListClass::Enumerate(ostream *out)
173{
174 periodentafel *periode = World::getInstance().getPeriode();
175 std::map<atomicNumber_t,unsigned int> counts;
176 double size=0;
177 Vector Origin;
178
179 // header
180 (*out) << "Index\tName\t\tAtoms\tFormula\tCenter\tSize" << endl;
181 (*out) << "-----------------------------------------------" << endl;
182 if (ListOfMolecules.size() == 0)
183 (*out) << "\tNone" << endl;
184 else {
185 Origin.Zero();
186 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
187 // count atoms per element and determine size of bounding sphere
188 size=0.;
189 for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
190 counts[(*iter)->getType()->getNumber()]++;
191 if ((*iter)->DistanceSquared(Origin) > size)
192 size = (*iter)->DistanceSquared(Origin);
193 }
194 // output Index, Name, number of atoms, chemical formula
195 (*out) << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->getAtomCount() << "\t";
196
197 std::map<atomicNumber_t,unsigned int>::reverse_iterator iter;
198 for(iter=counts.rbegin(); iter!=counts.rend();++iter){
199 atomicNumber_t Z =(*iter).first;
200 (*out) << periode->FindElement(Z)->getSymbol() << (*iter).second;
201 }
202 // Center and size
203 Vector *Center = (*ListRunner)->DetermineCenterOfAll();
204 (*out) << "\t" << *Center << "\t" << sqrt(size) << endl;
205 delete(Center);
206 }
207 }
208};
209
210/** Returns the molecule with the given index \a index.
211 * \param index index of the desired molecule
212 * \return pointer to molecule structure, NULL if not found
213 */
214molecule * MoleculeListClass::ReturnIndex(int index)
215{
216 for(MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
217 if ((*ListRunner)->IndexNr == index)
218 return (*ListRunner);
219 return NULL;
220};
221
222
223/** Simple output of the pointers in ListOfMolecules.
224 * \param *out output stream
225 */
226void MoleculeListClass::Output(ofstream *out)
227{
228 DoLog(1) && (Log() << Verbose(1) << "MoleculeList: ");
229 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
230 DoLog(0) && (Log() << Verbose(0) << *ListRunner << "\t");
231 DoLog(0) && (Log() << Verbose(0) << endl);
232};
233
234/** Calculates necessary hydrogen correction due to unwanted interaction between saturated ones.
235 * If for a pair of two hydrogen atoms a and b, at least is a saturated one, and a and b are not
236 * bonded to the same atom, then we add for this pair a correction term constructed from a Morse
237 * potential function fit to QM calculations with respecting to the interatomic hydrogen distance.
238 * \param &path path to file
239 */
240bool MoleculeListClass::AddHydrogenCorrection(std::string &path)
241{
242 bond *Binder = NULL;
243 double ***FitConstant = NULL, **correction = NULL;
244 int a, b;
245 ofstream output;
246 ifstream input;
247 string line;
248 stringstream zeile;
249 double distance;
250 char ParsedLine[1023];
251 double tmp;
252 char *FragmentNumber = NULL;
253
254 DoLog(1) && (Log() << Verbose(1) << "Saving hydrogen saturation correction ... ");
255 // 0. parse in fit constant files that should have the same dimension as the final energy files
256 // 0a. find dimension of matrices with constants
257 line = path;
258 line += "1";
259 line += FITCONSTANTSUFFIX;
260 input.open(line.c_str());
261 if (input.fail()) {
262 DoLog(1) && (Log() << Verbose(1) << endl << "Unable to open " << line << ", is the directory correct?" << endl);
263 return false;
264 }
265 a = 0;
266 b = -1; // we overcount by one
267 while (!input.eof()) {
268 input.getline(ParsedLine, 1023);
269 zeile.str(ParsedLine);
270 int i = 0;
271 while (!zeile.eof()) {
272 zeile >> distance;
273 i++;
274 }
275 if (i > a)
276 a = i;
277 b++;
278 }
279 DoLog(0) && (Log() << Verbose(0) << "I recognized " << a << " columns and " << b << " rows, ");
280 input.close();
281
282 // 0b. allocate memory for constants
283 FitConstant = new double**[3];
284 for (int k = 0; k < 3; k++) {
285 FitConstant[k] = new double*[a];
286 for (int i = a; i--;) {
287 FitConstant[k][i] = new double[b];
288 for (int j = b; j--;) {
289 FitConstant[k][i][j] = 0.;
290 }
291 }
292 }
293 // 0c. parse in constants
294 for (int i = 0; i < 3; i++) {
295 line = path;
296 line.append("/");
297 line += FRAGMENTPREFIX;
298 sprintf(ParsedLine, "%d", i + 1);
299 line += ParsedLine;
300 line += FITCONSTANTSUFFIX;
301 input.open(line.c_str());
302 if (input == NULL) {
303 DoeLog(0) && (eLog()<< Verbose(0) << endl << "Unable to open " << line << ", is the directory correct?" << endl);
304 performCriticalExit();
305 return false;
306 }
307 int k = 0, l;
308 while ((!input.eof()) && (k < b)) {
309 input.getline(ParsedLine, 1023);
310 //Log() << Verbose(0) << "Current Line: " << ParsedLine << endl;
311 zeile.str(ParsedLine);
312 zeile.clear();
313 l = 0;
314 while ((!zeile.eof()) && (l < a)) {
315 zeile >> FitConstant[i][l][k];
316 //Log() << Verbose(0) << FitConstant[i][l][k] << "\t";
317 l++;
318 }
319 //Log() << Verbose(0) << endl;
320 k++;
321 }
322 input.close();
323 }
324 for (int k = 0; k < 3; k++) {
325 DoLog(0) && (Log() << Verbose(0) << "Constants " << k << ":" << endl);
326 for (int j = 0; j < b; j++) {
327 for (int i = 0; i < a; i++) {
328 DoLog(0) && (Log() << Verbose(0) << FitConstant[k][i][j] << "\t");
329 }
330 DoLog(0) && (Log() << Verbose(0) << endl);
331 }
332 DoLog(0) && (Log() << Verbose(0) << endl);
333 }
334
335 // 0d. allocate final correction matrix
336 correction = new double*[a];
337 for (int i = a; i--;)
338 correction[i] = new double[b];
339
340 // 1a. go through every molecule in the list
341 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
342 // 1b. zero final correction matrix
343 for (int k = a; k--;)
344 for (int j = b; j--;)
345 correction[k][j] = 0.;
346 // 2. take every hydrogen that is a saturated one
347 for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
348 //Log() << Verbose(1) << "(*iter): " << *(*iter) << " with first bond " << *((*iter)->ListOfBonds.begin()) << "." << endl;
349 if (((*iter)->getType()->getAtomicNumber() == 1) && (((*iter)->father == NULL)
350 || ((*iter)->father->getType()->getAtomicNumber() != 1))) { // if it's a hydrogen
351 for (molecule::const_iterator runner = (*ListRunner)->begin(); runner != (*ListRunner)->end(); ++runner) {
352 //Log() << Verbose(2) << "Runner: " << *(*runner) << " with first bond " << *((*iter)->ListOfBonds.begin()) << "." << endl;
353 // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
354 Binder = *((*runner)->ListOfBonds.begin());
355 if (((*runner)->getType()->getAtomicNumber() == 1) && ((*runner)->nr > (*iter)->nr) && (Binder->GetOtherAtom((*runner)) != Binder->GetOtherAtom((*iter)))) { // (hydrogens have only one bonding partner!)
356 // 4. evaluate the morse potential for each matrix component and add up
357 distance = (*runner)->distance(*(*iter));
358 //Log() << Verbose(0) << "Fragment " << (*ListRunner)->name << ": " << *(*runner) << "<= " << distance << "=>" << *(*iter) << ":" << endl;
359 for (int k = 0; k < a; k++) {
360 for (int j = 0; j < b; j++) {
361 switch (k) {
362 case 1:
363 case 7:
364 case 11:
365 tmp = pow(FitConstant[0][k][j] * (1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]))), 2);
366 break;
367 default:
368 tmp = FitConstant[0][k][j] * pow(distance, FitConstant[1][k][j]) + FitConstant[2][k][j];
369 };
370 correction[k][j] -= tmp; // ground state is actually lower (disturbed by additional interaction)
371 //Log() << Verbose(0) << tmp << "\t";
372 }
373 //Log() << Verbose(0) << endl;
374 }
375 //Log() << Verbose(0) << endl;
376 }
377 }
378 }
379 }
380 // 5. write final matrix to file
381 line = path;
382 line.append("/");
383 line += FRAGMENTPREFIX;
384 FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), (*ListRunner)->IndexNr);
385 line += FragmentNumber;
386 delete[] (FragmentNumber);
387 line += HCORRECTIONSUFFIX;
388 output.open(line.c_str());
389 output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
390 for (int j = 0; j < b; j++) {
391 for (int i = 0; i < a; i++)
392 output << correction[i][j] << "\t";
393 output << endl;
394 }
395 output.close();
396 }
397 for (int i = a; i--;)
398 delete[](correction[i]);
399 delete[](correction);
400
401 line = path;
402 line.append("/");
403 line += HCORRECTIONSUFFIX;
404 output.open(line.c_str());
405 output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
406 for (int j = 0; j < b; j++) {
407 for (int i = 0; i < a; i++)
408 output << 0 << "\t";
409 output << endl;
410 }
411 output.close();
412 // 6. free memory of parsed matrices
413 for (int k = 0; k < 3; k++) {
414 for (int i = a; i--;) {
415 delete[](FitConstant[k][i]);
416 }
417 delete[](FitConstant[k]);
418 }
419 delete[](FitConstant);
420 DoLog(0) && (Log() << Verbose(0) << "done." << endl);
421 return true;
422};
423
424/** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
425 * \param &path path to file
426 * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
427 * \return true - file written successfully, false - writing failed
428 */
429bool MoleculeListClass::StoreForcesFile(std::string &path, int *SortIndex)
430{
431 bool status = true;
432 string filename(path);
433 filename += FORCESFILE;
434 ofstream ForcesFile(filename.c_str());
435 periodentafel *periode=World::getInstance().getPeriode();
436
437 // open file for the force factors
438 DoLog(1) && (Log() << Verbose(1) << "Saving force factors ... ");
439 if (!ForcesFile.fail()) {
440 //Log() << Verbose(1) << "Final AtomicForcesList: ";
441 //output << prefix << "Forces" << endl;
442 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
443 periodentafel::const_iterator elemIter;
444 for(elemIter=periode->begin();elemIter!=periode->end();++elemIter){
445 if ((*ListRunner)->hasElement((*elemIter).first)) { // if this element got atoms
446 for(molecule::iterator atomIter = (*ListRunner)->begin(); atomIter !=(*ListRunner)->end();++atomIter){
447 if ((*atomIter)->getType()->getNumber() == (*elemIter).first) {
448 if (((*atomIter)->GetTrueFather() != NULL) && ((*atomIter)->GetTrueFather() != (*atomIter))) {// if there is a rea
449 //Log() << Verbose(0) << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
450 ForcesFile << SortIndex[(*atomIter)->GetTrueFather()->nr] << "\t";
451 } else
452 // otherwise a -1 to indicate an added saturation hydrogen
453 ForcesFile << "-1\t";
454 }
455 }
456 }
457 }
458 ForcesFile << endl;
459 }
460 ForcesFile.close();
461 DoLog(1) && (Log() << Verbose(1) << "done." << endl);
462 } else {
463 status = false;
464 DoLog(1) && (Log() << Verbose(1) << "failed to open file " << filename << "." << endl);
465 }
466 ForcesFile.close();
467
468 return status;
469};
470
471/** Writes a config file for each molecule in the given \a **FragmentList.
472 * \param *out output stream for debugging
473 * \param &prefix path and prefix to the fragment config files
474 * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
475 * \return true - success (each file was written), false - something went wrong.
476 */
477bool MoleculeListClass::OutputConfigForListOfFragments(std::string &prefix, int *SortIndex)
478{
479 ofstream outputFragment;
480 std::string FragmentName;
481 char PathBackup[MAXSTRINGSIZE];
482 bool result = true;
483 bool intermediateResult = true;
484 Vector BoxDimension;
485 char *FragmentNumber = NULL;
486 char *path = NULL;
487 int FragmentCounter = 0;
488 ofstream output;
489 RealSpaceMatrix cell_size = World::getInstance().getDomain().getM();
490 RealSpaceMatrix cell_size_backup = cell_size;
491
492 // store the fragments as config and as xyz
493 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
494 // save default path as it is changed for each fragment
495 path = World::getInstance().getConfig()->GetDefaultPath();
496 if (path != NULL)
497 strcpy(PathBackup, path);
498 else {
499 DoeLog(0) && (eLog()<< Verbose(0) << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl);
500 performCriticalExit();
501 }
502
503 // correct periodic
504 (*ListRunner)->ScanForPeriodicCorrection();
505
506 // output xyz file
507 FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), FragmentCounter++);
508 FragmentName = prefix + FragmentNumber + ".conf.xyz";
509 outputFragment.open(FragmentName.c_str(), ios::out);
510 DoLog(2) && (Log() << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as XYZ ...");
511 if ((intermediateResult = (*ListRunner)->OutputXYZ(&outputFragment)))
512 DoLog(0) && (Log() << Verbose(0) << " done." << endl);
513 else
514 DoLog(0) && (Log() << Verbose(0) << " failed." << endl);
515 result = result && intermediateResult;
516 outputFragment.close();
517 outputFragment.clear();
518
519 // list atoms in fragment for debugging
520 DoLog(2) && (Log() << Verbose(2) << "Contained atoms: ");
521 for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
522 DoLog(0) && (Log() << Verbose(0) << (*iter)->getName() << " ");
523 }
524 DoLog(0) && (Log() << Verbose(0) << endl);
525
526 // center on edge
527 (*ListRunner)->CenterEdge(&BoxDimension);
528 (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
529 for (int k = 0; k < NDIM; k++) {
530 BoxDimension[k] = 2.5 * (World::getInstance().getConfig()->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
531 cell_size.at(k,k) = BoxDimension[k] * 2.;
532 }
533 World::getInstance().setDomain(cell_size);
534 (*ListRunner)->Translate(&BoxDimension);
535
536 // also calculate necessary orbitals
537 //(*ListRunner)->CalculateOrbitals(*World::getInstance().getConfig);
538
539 // change path in config
540 FragmentName = PathBackup;
541 FragmentName += "/";
542 FragmentName += FRAGMENTPREFIX;
543 FragmentName += FragmentNumber;
544 FragmentName += "/";
545 World::getInstance().getConfig()->SetDefaultPath(FragmentName.c_str());
546
547 // and save as config
548 FragmentName = prefix + FragmentNumber + ".conf";
549 DoLog(2) && (Log() << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as config ...");
550 if ((intermediateResult = World::getInstance().getConfig()->Save(FragmentName.c_str(), (*ListRunner)->elemente, (*ListRunner))))
551 DoLog(0) && (Log() << Verbose(0) << " done." << endl);
552 else
553 DoLog(0) && (Log() << Verbose(0) << " failed." << endl);
554 result = result && intermediateResult;
555
556 // restore old config
557 World::getInstance().getConfig()->SetDefaultPath(PathBackup);
558
559 // and save as mpqc input file
560 FragmentName = prefix + FragmentNumber + ".conf";
561 DoLog(2) && (Log() << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as mpqc input ...");
562 if ((intermediateResult = World::getInstance().getConfig()->SaveMPQC(FragmentName.c_str(), (*ListRunner))))
563 DoLog(2) && (Log() << Verbose(2) << " done." << endl);
564 else
565 DoLog(0) && (Log() << Verbose(0) << " failed." << endl);
566
567 result = result && intermediateResult;
568 //outputFragment.close();
569 //outputFragment.clear();
570 delete[](FragmentNumber);
571 }
572 DoLog(0) && (Log() << Verbose(0) << " done." << endl);
573
574 // printing final number
575 DoLog(2) && (Log() << Verbose(2) << "Final number of fragments: " << FragmentCounter << "." << endl);
576
577 // restore cell_size
578 World::getInstance().setDomain(cell_size_backup);
579
580 return result;
581};
582
583/** Counts the number of molecules with the molecule::ActiveFlag set.
584 * \return number of molecules with ActiveFlag set to true.
585 */
586int MoleculeListClass::NumberOfActiveMolecules()
587{
588 int count = 0;
589 for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
590 count += ((*ListRunner)->ActiveFlag ? 1 : 0);
591 return count;
592};
593
594/** Count all atoms in each molecule.
595 * \return number of atoms in the MoleculeListClass.
596 * TODO: the inner loop should be done by some (double molecule::CountAtom()) function
597 */
598int MoleculeListClass::CountAllAtoms() const
599{
600 int AtomNo = 0;
601 for (MoleculeList::const_iterator MolWalker = ListOfMolecules.begin(); MolWalker != ListOfMolecules.end(); MolWalker++) {
602 AtomNo += (*MolWalker)->size();
603 }
604 return AtomNo;
605}
606
607/***********
608 * Methods Moved here from the menus
609 */
610
611void MoleculeListClass::createNewMolecule(periodentafel *periode) {
612 OBSERVE;
613 molecule *mol = NULL;
614 mol = World::getInstance().createMolecule();
615 insert(mol);
616};
617
618void MoleculeListClass::loadFromXYZ(periodentafel *periode){
619 molecule *mol = NULL;
620 Vector center;
621 char filename[MAXSTRINGSIZE];
622 Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
623 mol = World::getInstance().createMolecule();
624 do {
625 Log() << Verbose(0) << "Enter file name: ";
626 cin >> filename;
627 } while (!mol->AddXYZFile(filename));
628 mol->SetNameFromFilename(filename);
629 // center at set box dimensions
630 mol->CenterEdge(&center);
631 RealSpaceMatrix domain;
632 for(int i =0;i<NDIM;++i)
633 domain.at(i,i) = center[i];
634 World::getInstance().setDomain(domain);
635 insert(mol);
636}
637
638void MoleculeListClass::setMoleculeFilename() {
639 char filename[MAXSTRINGSIZE];
640 int nr;
641 molecule *mol = NULL;
642 do {
643 Log() << Verbose(0) << "Enter index of molecule: ";
644 cin >> nr;
645 mol = ReturnIndex(nr);
646 } while (mol == NULL);
647 Log() << Verbose(0) << "Enter name: ";
648 cin >> filename;
649 mol->SetNameFromFilename(filename);
650}
651
652void MoleculeListClass::parseXYZIntoMolecule(){
653 char filename[MAXSTRINGSIZE];
654 int nr;
655 molecule *mol = NULL;
656 mol = NULL;
657 do {
658 Log() << Verbose(0) << "Enter index of molecule: ";
659 cin >> nr;
660 mol = ReturnIndex(nr);
661 } while (mol == NULL);
662 Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
663 do {
664 Log() << Verbose(0) << "Enter file name: ";
665 cin >> filename;
666 } while (!mol->AddXYZFile(filename));
667 mol->SetNameFromFilename(filename);
668};
669
670void MoleculeListClass::eraseMolecule(){
671 int nr;
672 molecule *mol = NULL;
673 Log() << Verbose(0) << "Enter index of molecule: ";
674 cin >> nr;
675 for(MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
676 if (nr == (*ListRunner)->IndexNr) {
677 mol = *ListRunner;
678 ListOfMolecules.erase(ListRunner);
679 World::getInstance().destroyMolecule(mol);
680 break;
681 }
682};
683
684
685/******************************************* Class MoleculeLeafClass ************************************************/
686
687/** Constructor for MoleculeLeafClass root leaf.
688 * \param *Up Leaf on upper level
689 * \param *PreviousLeaf NULL - We are the first leaf on this level, otherwise points to previous in list
690 */
691//MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *Up = NULL, MoleculeLeafClass *Previous = NULL)
692MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf = NULL) :
693 Leaf(NULL),
694 previous(PreviousLeaf)
695{
696 // if (Up != NULL)
697 // if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
698 // Up->DownLeaf = this;
699 // UpLeaf = Up;
700 // DownLeaf = NULL;
701 if (previous != NULL) {
702 MoleculeLeafClass *Walker = previous->next;
703 previous->next = this;
704 next = Walker;
705 } else {
706 next = NULL;
707 }
708};
709
710/** Destructor for MoleculeLeafClass.
711 */
712MoleculeLeafClass::~MoleculeLeafClass()
713{
714 // if (DownLeaf != NULL) {// drop leaves further down
715 // MoleculeLeafClass *Walker = DownLeaf;
716 // MoleculeLeafClass *Next;
717 // do {
718 // Next = Walker->NextLeaf;
719 // delete(Walker);
720 // Walker = Next;
721 // } while (Walker != NULL);
722 // // Last Walker sets DownLeaf automatically to NULL
723 // }
724 // remove the leaf itself
725 if (Leaf != NULL) {
726 World::getInstance().destroyMolecule(Leaf);
727 Leaf = NULL;
728 }
729 // remove this Leaf from level list
730 if (previous != NULL)
731 previous->next = next;
732 // } else { // we are first in list (connects to UpLeaf->DownLeaf)
733 // if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
734 // NextLeaf->UpLeaf = UpLeaf; // either null as we are top level or the upleaf of the first node
735 // if (UpLeaf != NULL)
736 // UpLeaf->DownLeaf = NextLeaf; // either null as we are only leaf or NextLeaf if we are just the first
737 // }
738 // UpLeaf = NULL;
739 if (next != NULL) // are we last in list
740 next->previous = previous;
741 next = NULL;
742 previous = NULL;
743};
744
745/** Adds \a molecule leaf to the tree.
746 * \param *ptr ptr to molecule to be added
747 * \param *Previous previous MoleculeLeafClass referencing level and which on the level
748 * \return true - success, false - something went wrong
749 */
750bool MoleculeLeafClass::AddLeaf(molecule *ptr, MoleculeLeafClass *Previous)
751{
752 return false;
753};
754
755/** Fills the bond structure of this chain list subgraphs that are derived from a complete \a *reference molecule.
756 * Calls this routine in each MoleculeLeafClass::next subgraph if it's not NULL.
757 * \param *out output stream for debugging
758 * \param *reference reference molecule with the bond structure to be copied
759 * \param **&ListOfLocalAtoms Lookup table for this subgraph and index of each atom in \a *reference, may be NULL on start, then it is filled
760 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
761 * \return true - success, false - faoilure
762 */
763bool MoleculeLeafClass::FillBondStructureFromReference(const molecule * const reference, atom **&ListOfLocalAtoms, bool FreeList)
764{
765 atom *OtherWalker = NULL;
766 atom *Father = NULL;
767 bool status = true;
768 int AtomNo;
769
770 DoLog(1) && (Log() << Verbose(1) << "Begin of FillBondStructureFromReference." << endl);
771 // fill ListOfLocalAtoms if NULL was given
772 if (!FillListOfLocalAtoms(ListOfLocalAtoms, reference->getAtomCount(), FreeList)) {
773 DoLog(1) && (Log() << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl);
774 return false;
775 }
776
777 if (status) {
778 DoLog(1) && (Log() << Verbose(1) << "Creating adjacency list for subgraph " << Leaf << "." << endl);
779 // remove every bond from the list
780 for(molecule::iterator AtomRunner = Leaf->begin(); AtomRunner != Leaf->end(); ++AtomRunner)
781 for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
782 if ((*BondRunner)->leftatom == *AtomRunner)
783 delete((*BondRunner));
784
785 for(molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end(); ++iter) {
786 Father = (*iter)->GetTrueFather();
787 AtomNo = Father->nr; // global id of the current walker
788 for (BondList::const_iterator Runner = Father->ListOfBonds.begin(); Runner != Father->ListOfBonds.end(); (++Runner)) {
789 OtherWalker = ListOfLocalAtoms[(*Runner)->GetOtherAtom((*iter)->GetTrueFather())->nr]; // local copy of current bond partner of walker
790 if (OtherWalker != NULL) {
791 if (OtherWalker->nr > (*iter)->nr)
792 Leaf->AddBond((*iter), OtherWalker, (*Runner)->BondDegree);
793 } else {
794 DoLog(1) && (Log() << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << (*Runner)->GetOtherAtom((*iter)->GetTrueFather())->nr << "] is NULL!" << endl);
795 status = false;
796 }
797 }
798 }
799 }
800
801 if ((FreeList) && (ListOfLocalAtoms != NULL)) {
802 // free the index lookup list
803 delete[](ListOfLocalAtoms);
804 }
805 DoLog(1) && (Log() << Verbose(1) << "End of FillBondStructureFromReference." << endl);
806 return status;
807};
808
809/** Fills the root stack for sites to be used as root in fragmentation depending on order or adaptivity criteria
810 * Again, as in \sa FillBondStructureFromReference steps recursively through each Leaf in this chain list of molecule's.
811 * \param *out output stream for debugging
812 * \param *&RootStack stack to be filled
813 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site
814 * \param &FragmentCounter counts through the fragments in this MoleculeLeafClass
815 * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
816 */
817bool MoleculeLeafClass::FillRootStackForSubgraphs(KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
818{
819 atom *Father = NULL;
820
821 if (RootStack != NULL) {
822 // find first root candidates
823 if (&(RootStack[FragmentCounter]) != NULL) {
824 RootStack[FragmentCounter].clear();
825 for(molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end(); ++iter) {
826 Father = (*iter)->GetTrueFather();
827 if (AtomMask[Father->nr]) // apply mask
828#ifdef ADDHYDROGEN
829 if ((*iter)->getType()->getAtomicNumber() != 1) // skip hydrogen
830#endif
831 RootStack[FragmentCounter].push_front((*iter)->nr);
832 }
833 if (next != NULL)
834 next->FillRootStackForSubgraphs(RootStack, AtomMask, ++FragmentCounter);
835 } else {
836 DoLog(1) && (Log() << Verbose(1) << "Rootstack[" << FragmentCounter << "] is NULL." << endl);
837 return false;
838 }
839 FragmentCounter--;
840 return true;
841 } else {
842 DoLog(1) && (Log() << Verbose(1) << "Rootstack is NULL." << endl);
843 return false;
844 }
845};
846
847/** Fills a lookup list of father's Atom::nr -> atom for each subgraph.
848 * \param *out output stream from debugging
849 * \param **&ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
850 * \param GlobalAtomCount number of atoms in the complete molecule
851 * \param &FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
852 * \return true - success, false - failure (ListOfLocalAtoms != NULL)
853 */
854bool MoleculeLeafClass::FillListOfLocalAtoms(atom **&ListOfLocalAtoms, const int GlobalAtomCount, bool &FreeList)
855{
856 bool status = true;
857
858 if (ListOfLocalAtoms == NULL) { // allocate and fill list of this fragment/subgraph
859 status = status && Leaf->CreateFatherLookupTable(ListOfLocalAtoms, GlobalAtomCount);
860 FreeList = FreeList && true;
861 } else
862 return false;
863
864 return status;
865};
866
867/** The indices per keyset are compared to the respective father's Atom::nr in each subgraph and thus put into \a **&FragmentList.
868 * \param *out output stream fro debugging
869 * \param *reference reference molecule with the bond structure to be copied
870 * \param *KeySetList list with all keysets
871 * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
872 * \param **&FragmentList list to be allocated and returned
873 * \param &FragmentCounter counts the fragments as we move along the list
874 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
875 * \retuen true - success, false - failure
876 */
877bool MoleculeLeafClass::AssignKeySetsToFragment(molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList)
878{
879 bool status = true;
880 int KeySetCounter = 0;
881
882 DoLog(1) && (Log() << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl);
883 // fill ListOfLocalAtoms if NULL was given
884 if (!FillListOfLocalAtoms(ListOfLocalAtoms[FragmentCounter], reference->getAtomCount(), FreeList)) {
885 DoLog(1) && (Log() << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl);
886 return false;
887 }
888
889 // allocate fragment list
890 if (FragmentList == NULL) {
891 KeySetCounter = Count();
892 FragmentList = new Graph*[KeySetCounter];
893 for (int i=0;i<KeySetCounter;i++)
894 FragmentList[i] = NULL;
895 KeySetCounter = 0;
896 }
897
898 if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
899 // assign scanned keysets
900 if (FragmentList[FragmentCounter] == NULL)
901 FragmentList[FragmentCounter] = new Graph;
902 KeySet *TempSet = new KeySet;
903 for (Graph::iterator runner = KeySetList->begin(); runner != KeySetList->end(); runner++) { // key sets contain global numbers!
904 if (ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
905 // translate keyset to local numbers
906 for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
907 TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
908 // insert into FragmentList
909 FragmentList[FragmentCounter]->insert(GraphPair(*TempSet, pair<int, double> (KeySetCounter++, (*runner).second.second)));
910 }
911 TempSet->clear();
912 }
913 delete (TempSet);
914 if (KeySetCounter == 0) {// if there are no keysets, delete the list
915 DoLog(1) && (Log() << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl);
916 delete (FragmentList[FragmentCounter]);
917 } else
918 DoLog(1) && (Log() << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl);
919 FragmentCounter++;
920 if (next != NULL)
921 next->AssignKeySetsToFragment(reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
922 FragmentCounter--;
923 } else
924 DoLog(1) && (Log() << Verbose(1) << "KeySetList is NULL or empty." << endl);
925
926 if ((FreeList) && (ListOfLocalAtoms != NULL)) {
927 // free the index lookup list
928 delete[](ListOfLocalAtoms[FragmentCounter]);
929 }
930 DoLog(1) && (Log() << Verbose(1) << "End of AssignKeySetsToFragment." << endl);
931 return status;
932};
933
934/** Translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
935 * \param *out output stream for debugging
936 * \param **FragmentList Graph with local numbers per fragment
937 * \param &FragmentCounter counts the fragments as we move along the list
938 * \param &TotalNumberOfKeySets global key set counter
939 * \param &TotalGraph Graph to be filled with global numbers
940 */
941void MoleculeLeafClass::TranslateIndicesToGlobalIDs(Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph)
942{
943 DoLog(1) && (Log() << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl);
944 KeySet *TempSet = new KeySet;
945 if (FragmentList[FragmentCounter] != NULL) {
946 for (Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
947 for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
948 TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
949 TotalGraph.insert(GraphPair(*TempSet, pair<int, double> (TotalNumberOfKeySets++, (*runner).second.second)));
950 TempSet->clear();
951 }
952 delete (TempSet);
953 } else {
954 DoLog(1) && (Log() << Verbose(1) << "FragmentList is NULL." << endl);
955 }
956 if (next != NULL)
957 next->TranslateIndicesToGlobalIDs(FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
958 FragmentCounter--;
959 DoLog(1) && (Log() << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl);
960};
961
962/** Simply counts the number of items in the list, from given MoleculeLeafClass.
963 * \return number of items
964 */
965int MoleculeLeafClass::Count() const
966{
967 if (next != NULL)
968 return next->Count() + 1;
969 else
970 return 1;
971};
972
Note: See TracBrowser for help on using the repository browser.