source: src/periodentafel.cpp@ 2f40c0e

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 2f40c0e was 274d45, checked in by Frederik Heber <heber@…>, 15 years ago

FIX: Atoms were stored not in the sequence they were loaded.

  1. The main problem is molecule::atomSet which is a set<atom *>, i.e. atoms are sorted by their appearance in memory. As memory need not be allocated sequentially, this gives rise to extreme arbitririty which is not desired. Instead the atoms should be stored in the sequence they were loaded/created. The solution is as follows:
  • config::SaveAll()
  • molecule::atomSet is now a list<atom *>
  • molecule::atomIds is a new set<atomId_t> (atomIdSet) which controls that (global) ids remain unique in the no more Atomset's set (but list)
  • molecule::erase() erases also in atomIds
  • molecule::insert() checks whether id is present by atomIds
  • molecule::find() as std::list does not have a find, we just go through the list until the object is found (or not), this may be speeded up by another internal list.
  • molecule::InternalPointer made lots of confusion as virtual function GoToFirst() is const, hence begin() (needed therein) returns const_iterator, which then cannot be simply re-cast into an iterator: We make it a pointer, reinterpret_cast the pointer and reference it back. Although InternalPointer is mutable, the compiler cannot use the non-const function begin() (it cannot be made const, as overloading is not allowed). (this is noted in the code also extensively.)
  • molecule::containsAtom() does not use count but the new find, as it returns only boolean anyway.
  • rewrote MoleculeListClass::SimpleMerge() to get rid of the extra iterator, as we remove all atoms in the end anyway.
  • FIX: MoleculeListClass::SimpleMultiMerge() - the created mol was not inserted into the moleculelist in the end, although it is the only one to remain.
  1. All other databases had missing headers with respect to those stored in elements_db.cpp. Hence, valence of hydrogen was not parsed and this caused several failures in CalculateOrbitals() (Psi numbers and MaxMinSetp in pcp conf file).
  1. Subsequenytly, various test cases (12, 21, 30, 31, 36-38, 39) were broken. This had two reasons:
  • Seemingly, CalculateOrbitals() was broken before hence MaxMinStep (pcp config) and MaxPsiDouble/PsiMaxNn[Up|Down] were always 0. (10-21,30-31,39)
  • As the order is now correct, fixes from commits c9217161ec2a5d5db508557fe98a32068461f45b and 22a6da8380911571debebd69444d2615450bbbd8 were obselete and have been reverted (order of the Ion?_Type...): Molecules/6 (30), Molecules/7 (31), Filling/1 (39)
  • Due to different ordering, Tesselation/3 (38) had completely different .dat file (though same tesselation)
  • r3d had small differences, mostly order or epsilon (0 not being 0 but ..-e16), hence diffing was deactivated, as r3d is deprecated anyway (since vmd can render triangles as well and better).

Signed-off-by: Frederik Heber <heber@…>

  • Property mode set to 100755
File size: 15.1 KB
Line 
1/** \file periodentafel.cpp
2 *
3 * Function implementations for the class periodentafel.
4 *
5 */
6
7#include "Helpers/MemDebug.hpp"
8
9using namespace std;
10
11#include <iomanip>
12#include <iostream>
13#include <fstream>
14#include <cstring>
15
16#include "Helpers/Assert.hpp"
17#include "element.hpp"
18#include "elements_db.hpp"
19#include "helpers.hpp"
20#include "lists.hpp"
21#include "log.hpp"
22#include "periodentafel.hpp"
23#include "verbose.hpp"
24
25using namespace std;
26
27/************************************* Functions for class periodentafel ***************************/
28
29/** constructor for class periodentafel
30 * Initialises start and end of list and resets periodentafel::checkliste to false.
31 */
32periodentafel::periodentafel()
33{
34 bool status = true;
35 status = LoadElementsDatabase(new stringstream(elementsDB,ios_base::in));
36 ASSERT(status, "General element initialization failed");
37 status = LoadValenceDatabase(new stringstream(valenceDB,ios_base::in));
38 ASSERT(status, "Valence entry of element initialization failed");
39 status = LoadOrbitalsDatabase(new stringstream(orbitalsDB,ios_base::in));
40 ASSERT(status, "Orbitals entry of element initialization failed");
41 status = LoadHBondAngleDatabase(new stringstream(HbondangleDB,ios_base::in));
42 ASSERT(status, "HBond angle entry of element initialization failed");
43 status = LoadHBondLengthsDatabase(new stringstream(HbonddistanceDB,ios_base::in));
44 ASSERT(status, "HBond distance entry of element initialization failed");
45};
46
47/** destructor for class periodentafel
48 * Removes every element and afterwards deletes start and end of list.
49 * TODO: Handle when elements have changed and store databases then
50 */
51periodentafel::~periodentafel()
52{
53 CleanupPeriodtable();
54};
55
56/** Adds element to period table list
57 * \param *pointer element to be added
58 * \return iterator to added element
59 */
60periodentafel::iterator periodentafel::AddElement(element * const pointer)
61{
62 atomicNumber_t Z = pointer->getNumber();
63 ASSERT(!elements.count(Z), "Element is already present.");
64 pointer->sort = &pointer->Z;
65 if (pointer->getNumber() < 1 && pointer->getNumber() >= MAX_ELEMENTS)
66 DoeLog(0) && (eLog() << Verbose(0) << "Invalid Z number!\n");
67 pair<iterator,bool> res = elements.insert(pair<atomicNumber_t,element*>(Z,pointer));
68 return res.first;
69};
70
71/** Removes element from list.
72 * \param *pointer element to be removed
73 */
74size_t periodentafel::RemoveElement(element * const pointer)
75{
76 return RemoveElement(pointer->getNumber());
77};
78
79/** Removes element from list.
80 * \param Z element to be removed
81 */
82size_t periodentafel::RemoveElement(atomicNumber_t Z)
83{
84 return elements.erase(Z);
85};
86
87/** Removes every element from the period table.
88 */
89void periodentafel::CleanupPeriodtable()
90{
91 for(iterator iter=elements.begin();iter!=elements.end();++iter){
92 delete(*iter).second;
93 }
94 elements.clear();
95};
96
97/** Finds an element by its atomic number.
98 * If element is not yet in list, returns NULL.
99 * \param Z atomic number
100 * \return pointer to element or NULL if not found
101 */
102element * const periodentafel::FindElement(atomicNumber_t Z) const
103{
104 const_iterator res = elements.find(Z);
105 return res!=elements.end()?((*res).second):0;
106};
107
108/** Finds an element by its atomic number.
109 * If element is not yet in list, datas are asked and stored in database.
110 * \param shorthand chemical symbol of the element, e.g. H for hydrogene
111 * \return pointer to element
112 */
113element * const periodentafel::FindElement(const char * const shorthand) const
114{
115 element *res = 0;
116 for(const_iterator iter=elements.begin();iter!=elements.end();++iter) {
117 if((*iter).second->getSymbol() == shorthand){
118 res = (*iter).second;
119 break;
120 }
121 }
122 return res;
123};
124
125/** Asks for element number and returns pointer to element
126 * \return desired element or NULL
127 */
128element * const periodentafel::AskElement() const
129{
130 element * walker = NULL;
131 int Z;
132 do {
133 DoLog(0) && (Log() << Verbose(0) << "Atomic number Z: ");
134 cin >> Z;
135 walker = this->FindElement(Z); // give type
136 } while (walker == NULL);
137 return walker;
138};
139
140/** Asks for element and if not found, presents mask to enter info.
141 * \return pointer to either present or newly created element
142 */
143element * const periodentafel::EnterElement()
144{
145 atomicNumber_t Z = 0;
146 DoLog(0) && (Log() << Verbose(0) << "Atomic number: " << Z << endl);
147 cin >> Z;
148 element * const res = FindElement(Z);
149 if (!res) {
150 // TODO: make this using the constructor
151 DoLog(0) && (Log() << Verbose(0) << "Element not found in database, please enter." << endl);
152 element *tmp = new element;
153 tmp->Z = Z;
154 DoLog(0) && (Log() << Verbose(0) << "Mass: " << endl);
155 cin >> tmp->mass;
156 DoLog(0) && (Log() << Verbose(0) << "Name [max 64 chars]: " << endl);
157 cin >> tmp->name;
158 DoLog(0) && (Log() << Verbose(0) << "Short form [max 3 chars]: " << endl);
159 cin >> tmp->symbol;
160 AddElement(tmp);
161 return tmp;
162 }
163 return res;
164};
165
166
167/******************** Access to iterators ****************************/
168periodentafel::const_iterator periodentafel::begin(){
169 return elements.begin();
170}
171
172periodentafel::const_iterator periodentafel::end(){
173 return elements.end();
174}
175
176periodentafel::reverse_iterator periodentafel::rbegin(){
177 return reverse_iterator(elements.end());
178}
179
180periodentafel::reverse_iterator periodentafel::rend(){
181 return reverse_iterator(elements.begin());
182}
183
184/** Prints period table to given stream.
185 * \param output stream
186 */
187bool periodentafel::Output(ostream * const output) const
188{
189 bool result = true;
190 if (output != NULL) {
191 for(const_iterator iter=elements.begin(); iter !=elements.end();++iter){
192 result = result && (*iter).second->Output(output);
193 }
194 return result;
195 } else
196 return false;
197};
198
199/** Prints period table to given stream.
200 * \param *output output stream
201 * \param *checkliste elements table for this molecule
202 */
203bool periodentafel::Checkout(ostream * const output, const int * const checkliste) const
204{
205 bool result = true;
206 int No = 1;
207
208 if (output != NULL) {
209 *output << "# Ion type data (PP = PseudoPotential, Z = atomic number)" << endl;
210 *output << "#Ion_TypeNr.\tAmount\tZ\tRGauss\tL_Max(PP)L_Loc(PP)IonMass\t# chemical name, symbol" << endl;
211 for(const_iterator iter=elements.begin(); iter!=elements.end();++iter){
212 if (((*iter).first < MAX_ELEMENTS) && (checkliste[(*iter).first])) {
213 (*iter).second->No = No;
214 result = result && (*iter).second->Checkout(output, No++, checkliste[(*iter).first]);
215 }
216 }
217 return result;
218 } else
219 return false;
220};
221
222/** Loads element list from file.
223 * \param *path to to standard file names
224 */
225bool periodentafel::LoadPeriodentafel(const char *path)
226{
227 ifstream input;
228 bool status = true;
229 bool otherstatus = true;
230 char filename[255];
231
232 // fill elements DB
233 strncpy(filename, path, MAXSTRINGSIZE);
234 strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
235 strncat(filename, STANDARDELEMENTSDB, MAXSTRINGSIZE-strlen(filename));
236 input.open(filename);
237 if (!input.fail())
238 DoLog(0) && (Log() << Verbose(0) << "Using " << filename << " as elements database." << endl);
239 status = status && LoadElementsDatabase(&input);
240 input.close();
241 input.clear();
242
243 // fill valence DB per element
244 strncpy(filename, path, MAXSTRINGSIZE);
245 strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
246 strncat(filename, STANDARDVALENCEDB, MAXSTRINGSIZE-strlen(filename));
247 input.open(filename);
248 if (!input.fail())
249 DoLog(0) && (Log() << Verbose(0) << "Using " << filename << " as valence database." << endl);
250 otherstatus = otherstatus && LoadValenceDatabase(&input);
251 input.close();
252 input.clear();
253
254 // fill orbitals DB per element
255 strncpy(filename, path, MAXSTRINGSIZE);
256 strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
257 strncat(filename, STANDARDORBITALDB, MAXSTRINGSIZE-strlen(filename));
258 input.open(filename);
259 if (!input.fail())
260 DoLog(0) && (Log() << Verbose(0) << "Using " << filename << " as orbitals database." << endl);
261 otherstatus = otherstatus && LoadOrbitalsDatabase(&input);
262 input.close();
263 input.clear();
264
265 // fill H-BondAngle DB per element
266 strncpy(filename, path, MAXSTRINGSIZE);
267 strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
268 strncat(filename, STANDARDHBONDANGLEDB, MAXSTRINGSIZE-strlen(filename));
269 input.open(filename);
270 if (!input.fail())
271 DoLog(0) && (Log() << Verbose(0) << "Using " << filename << " as H bond angle database." << endl);
272 otherstatus = otherstatus && LoadHBondAngleDatabase(&input);
273 input.close();
274 input.clear();
275
276 // fill H-BondDistance DB per element
277 strncpy(filename, path, MAXSTRINGSIZE);
278 strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
279 strncat(filename, STANDARDHBONDDISTANCEDB, MAXSTRINGSIZE-strlen(filename));
280 input.open(filename);
281 if (!input.fail())
282 DoLog(0) && (Log() << Verbose(0) << "Using " << filename << " as H bond length database." << endl);
283 otherstatus = otherstatus && LoadHBondLengthsDatabase(&input);
284 input.close();
285 input.clear();
286
287 if (!otherstatus){
288 DoeLog(2) && (eLog()<< Verbose(2) << "Something went wrong while parsing the other databases!" << endl);
289 }
290
291 return status;
292};
293
294/** load the element info.
295 * \param *input stream to parse from
296 * \return true - parsing successful, false - something went wrong
297 */
298bool periodentafel::LoadElementsDatabase(istream *input)
299{
300 bool status = true;
301 int counter = 0;
302 pair< std::map<atomicNumber_t,element*>::iterator, bool > InserterTest;
303 if (!(*input).fail()) {
304 (*input).getline(header1, MAXSTRINGSIZE);
305 (*input).getline(header2, MAXSTRINGSIZE); // skip first two header lines
306 DoLog(0) && (Log() << Verbose(0) << "Parsed elements:");
307 while (!(*input).eof()) {
308 element *neues = new element;
309 (*input) >> neues->name;
310 //(*input) >> ws;
311 (*input) >> neues->symbol;
312 //(*input) >> ws;
313 (*input) >> neues->period;
314 //(*input) >> ws;
315 (*input) >> neues->group;
316 //(*input) >> ws;
317 (*input) >> neues->block;
318 //(*input) >> ws;
319 (*input) >> neues->Z;
320 //(*input) >> ws;
321 (*input) >> neues->mass;
322 //(*input) >> ws;
323 (*input) >> neues->CovalentRadius;
324 //(*input) >> ws;
325 (*input) >> neues->VanDerWaalsRadius;
326 //(*input) >> ws;
327 (*input) >> ws;
328 //neues->Output((ofstream *)&cout);
329 if ((neues->getNumber() > 0) && (neues->getNumber() < MAX_ELEMENTS)) {
330 if (elements.count(neues->getNumber())) {// if element already present, remove and delete old one (i.e. replace it)
331 //cout << neues->symbol << " is present already." << endl;
332 element * const Elemental = FindElement(neues->getNumber());
333 ASSERT(Elemental != NULL, "element should be present but is not??");
334 *Elemental = *neues;
335 } else {
336 InserterTest = elements.insert(pair <atomicNumber_t,element*> (neues->getNumber(), neues));
337 ASSERT(InserterTest.second, "Could not insert new element into periodentafel on LoadElementsDatabase().");
338 }
339 DoLog(0) && (Log() << Verbose(0) << " " << elements[neues->getNumber()]->symbol);
340 counter++;
341 } else {
342 DoeLog(2) && (eLog() << Verbose(2) << "Detected empty line or invalid element in elements db, discarding." << endl);
343 DoLog(0) && (Log() << Verbose(0) << " <?>");
344 delete(neues);
345 }
346 }
347 DoLog(0) && (Log() << Verbose(0) << endl);
348 } else {
349 DoeLog(1) && (eLog() << Verbose(1) << "Could not open the database." << endl);
350 status = false;
351 }
352
353 if (counter == 0)
354 status = false;
355
356 return status;
357}
358
359/** load the valence info.
360 * \param *input stream to parse from
361 * \return true - parsing successful, false - something went wrong
362 */
363bool periodentafel::LoadValenceDatabase(istream *input)
364{
365 char dummy[MAXSTRINGSIZE];
366 if (!(*input).fail()) {
367 (*input).getline(dummy, MAXSTRINGSIZE);
368 while (!(*input).eof()) {
369 atomicNumber_t Z;
370 (*input) >> Z;
371 ASSERT(elements.count(Z), "Element not present");
372 (*input) >> ws;
373 (*input) >> elements[Z]->Valence;
374 (*input) >> ws;
375 //Log() << Verbose(3) << "Element " << Z << " has " << FindElement(Z)->Valence << " valence electrons." << endl;
376 }
377 return true;
378 } else
379 return false;
380}
381
382/** load the orbitals info.
383 * \param *input stream to parse from
384 * \return true - parsing successful, false - something went wrong
385 */
386bool periodentafel::LoadOrbitalsDatabase(istream *input)
387{
388 char dummy[MAXSTRINGSIZE];
389 if (!(*input).fail()) {
390 (*input).getline(dummy, MAXSTRINGSIZE);
391 while (!(*input).eof()) {
392 atomicNumber_t Z;
393 (*input) >> Z;
394 ASSERT(elements.count(Z), "Element not present");
395 (*input) >> ws;
396 (*input) >> elements[Z]->NoValenceOrbitals;
397 (*input) >> ws;
398 //Log() << Verbose(3) << "Element " << Z << " has " << FindElement(Z)->NoValenceOrbitals << " number of singly occupied valence orbitals." << endl;
399 }
400 return true;
401 } else
402 return false;
403}
404
405/** load the hbond angles info.
406 * \param *input stream to parse from
407 * \return true - parsing successful, false - something went wrong
408 */
409bool periodentafel::LoadHBondAngleDatabase(istream *input)
410{
411 char dummy[MAXSTRINGSIZE];
412 if (!(*input).fail()) {
413 (*input).getline(dummy, MAXSTRINGSIZE);
414 while (!(*input).eof()) {
415 atomicNumber_t Z;
416 (*input) >> Z;
417 ASSERT(elements.count(Z), "Element not present");
418 (*input) >> ws;
419 (*input) >> elements[Z]->HBondAngle[0];
420 (*input) >> elements[Z]->HBondAngle[1];
421 (*input) >> elements[Z]->HBondAngle[2];
422 (*input) >> ws;
423 //Log() << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->HBondAngle[0] << ", " << FindElement((int)tmp)->HBondAngle[1] << ", " << FindElement((int)tmp)->HBondAngle[2] << " degrees bond angle for one, two, three connected hydrogens." << endl;
424 }
425 return true;
426 } else
427 return false;
428}
429
430/** load the hbond lengths info.
431 * \param *input stream to parse from
432 * \return true - parsing successful, false - something went wrong
433 */
434bool periodentafel::LoadHBondLengthsDatabase(istream *input)
435{
436 char dummy[MAXSTRINGSIZE];
437 if (!(*input).fail()) {
438 (*input).getline(dummy, MAXSTRINGSIZE);
439 while (!(*input).eof()) {
440 atomicNumber_t Z;
441 (*input) >> Z;
442 ASSERT(elements.count(Z), "Element not present");
443 (*input) >> ws;
444 (*input) >> elements[Z]->HBondDistance[0];
445 (*input) >> elements[Z]->HBondDistance[1];
446 (*input) >> elements[Z]->HBondDistance[2];
447 (*input) >> ws;
448 //Log() << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->HBondDistance[0] << " Angstrom typical distance to hydrogen." << endl;
449 }
450 return true;
451 } else
452 return false;
453}
454
455/** Stores element list to file.
456 */
457bool periodentafel::StorePeriodentafel(const char *path) const
458{
459 bool result = true;
460 ofstream f;
461 char filename[MAXSTRINGSIZE];
462
463 strncpy(filename, path, MAXSTRINGSIZE);
464 strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
465 strncat(filename, STANDARDELEMENTSDB, MAXSTRINGSIZE-strlen(filename));
466 f.open(filename);
467 if (f != NULL) {
468 f << header1 << endl;
469 f << header2 << endl;
470 for(const_iterator iter=elements.begin();iter!=elements.end();++iter){
471 result = result && (*iter).second->Output(&f);
472 }
473 f.close();
474 return true;
475 } else
476 return result;
477};
Note: See TracBrowser for help on using the repository browser.