source: src/Formula.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 94d5ac6, checked in by Frederik Heber <heber@…>, 12 years ago

FIX: As we use GSL internally, we are as of now required to use GPL v2 license.

  • GNU Scientific Library is used at every place in the code, especially the sub-package LinearAlgebra is based on it which in turn is used really everywhere in the remainder of MoleCuilder. Hence, we have to use the GPL license for the whole of MoleCuilder. In effect, GPL's COPYING was present all along and stated the terms of the GPL v2 license.
  • Hence, I added the default GPL v2 disclaimer to every source file and removed the note about a (actually missing) LICENSE file.
  • also, I added a help-redistribute action which again gives the disclaimer of the GPL v2.
  • also, I changed in the disclaimer that is printed at every program start in builder_init.cpp.
  • TEST: Added check on GPL statement present in every module to test CodeChecks project-disclaimer.
  • Property mode set to 100644
File size: 15.1 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/*
24 * Formula.cpp
25 *
26 * Created on: Jul 21, 2010
27 * Author: crueger
28 */
29
30// include config.h
31#ifdef HAVE_CONFIG_H
32#include <config.h>
33#endif
34
35#include "CodePatterns/MemDebug.hpp"
36
37#include "Formula.hpp"
38
39#include <sstream>
40
41#include "World.hpp"
42#include "Element/periodentafel.hpp"
43#include "Element/element.hpp"
44#include "CodePatterns/Assert.hpp"
45#include "CodePatterns/Range.hpp"
46
47using namespace std;
48
49Formula::Formula() :
50 numElements(0)
51{}
52
53Formula::Formula(const Formula &src) :
54 elementCounts(src.elementCounts),
55 numElements(src.numElements)
56{}
57
58Formula::Formula(const string &formula) :
59 numElements(0)
60{
61 fromString(formula);
62}
63
64Formula::~Formula()
65{}
66
67Formula &Formula::operator=(const Formula &rhs){
68 // No self-assignment check needed
69 elementCounts=rhs.elementCounts;
70 numElements=rhs.numElements;
71 return *this;
72}
73
74std::string Formula::toString() const{
75 stringstream sstr;
76 for(const_iterator iter=end();iter!=begin();){
77 --iter;
78 sstr << (*iter).first->getSymbol();
79 if((*iter).second>1)
80 sstr << (*iter).second;
81 }
82 return sstr.str();
83}
84
85void Formula::fromString(const std::string &formula) throw(FormulaStringParseException){
86 // make this transactional, in case an error is thrown
87 Formula res;
88 string::const_iterator begin = formula.begin();
89 string::const_iterator end = formula.end();
90 res.parseFromString(begin,end,static_cast<char>(0));
91 (*this)=res;
92}
93
94int Formula::parseMaybeNumber(std::string::const_iterator &it,string::const_iterator &end) throw(FormulaStringParseException){
95 static const range<char> Numbers = makeRange('0',static_cast<char>('9'+1));
96 int count = 0;
97 while(it!=end && Numbers.isInRange(*it))
98 count = (count*10) + ((*it++)-Numbers.first);
99 // one is implicit
100 count = (count!=0)?count:1;
101 return count;
102}
103
104void Formula::parseFromString(std::string::const_iterator &it,string::const_iterator &end,char delimiter) throw(FormulaStringParseException){
105 // some constants needed for parsing... Assumes ASCII, change if other encodings are used
106 static const range<char> CapitalLetters = makeRange('A',static_cast<char>('Z'+1));
107 static const range<char> SmallLetters = makeRange('a',static_cast<char>('z'+1));
108 map<char,char> delimiters;
109 delimiters['('] = ')';
110 delimiters['['] = ']';
111 // clean the formula
112 clear();
113 for(/*send from above*/;it!=end && *it!=delimiter;/*updated in loop*/){
114 // we might have a sub formula
115 if(delimiters.count(*it)){
116 Formula sub;
117 char nextdelim=delimiters[*it];
118 sub.parseFromString(++it,end,nextdelim);
119 if(!sub.getElementCount()){
120 throw FormulaStringParseException() << FormulaString( string(it, end) );
121 }
122 int count = parseMaybeNumber(++it,end);
123 addFormula(sub,count);
124 continue;
125 }
126 string shorthand;
127 // Atom names start with a capital letter
128 if(!CapitalLetters.isInRange(*it))
129 throw FormulaStringParseException() << FormulaString( string(it, end) );
130 shorthand+=(*it++);
131 // the rest of the name follows
132 while(it!=end && SmallLetters.isInRange(*it))
133 shorthand+=(*it++);
134 int count = parseMaybeNumber(it,end);
135 // test if the shorthand exists
136 if(!World::getInstance().getPeriode()->FindElement(shorthand))
137 throw FormulaStringParseException() << FormulaString( string(it, end) );
138 // done, we can get the next one
139 addElements(shorthand,count);
140 }
141 if(it==end && delimiter!=0){
142 throw FormulaStringParseException() << FormulaString( string(it, end) );
143 }
144}
145
146unsigned int Formula::getElementCount() const{
147 return numElements;
148}
149
150bool Formula::hasElement(const element *element) const{
151 ASSERT(element,"Invalid pointer in Formula::hasElement(element*)");
152 return hasElement(element->getAtomicNumber());
153}
154
155bool Formula::hasElement(atomicNumber_t Z) const{
156 ASSERT(Z>0,"Invalid atomic Number");
157 ASSERT(World::getInstance().getPeriode()->FindElement(Z),"No Element with this number in Periodentafel");
158 return elementCounts.size()>=Z && elementCounts[Z-1];
159}
160
161bool Formula::hasElement(const string &shorthand) const{
162 const element * element = World::getInstance().getPeriode()->FindElement(shorthand);
163 return hasElement(element);
164}
165
166void Formula::operator+=(const element *element){
167 ASSERT(element,"Invalid pointer in increment of Formula");
168 operator+=(element->getAtomicNumber());
169}
170
171void Formula::operator+=(atomicNumber_t Z){
172 ASSERT(Z>0,"Invalid atomic Number");
173 ASSERT(World::getInstance().getPeriode()->FindElement(Z),"No Element with this number in Periodentafel");
174 elementCounts.resize(max<atomicNumber_t>(Z,elementCounts.size()),0); // No-op when we already have the right size
175 // might need to update number of elements
176 if(!elementCounts[Z-1]){
177 numElements++;
178 }
179 elementCounts[Z-1]++; // atomic numbers start at 1
180}
181
182void Formula::operator+=(const string &shorthand){
183 const element * element = World::getInstance().getPeriode()->FindElement(shorthand);
184 operator+=(element);
185}
186
187void Formula::operator-=(const element *element){
188 ASSERT(element,"Invalid pointer in decrement of Formula");
189 operator-=(element->getAtomicNumber());
190}
191
192void Formula::operator-=(atomicNumber_t Z){
193 ASSERT(Z>0,"Invalid atomic Number");
194 ASSERT(World::getInstance().getPeriode()->FindElement(Z),"No Element with this number in Periodentafel");
195 ASSERT(elementCounts.size()>=Z && elementCounts[Z-1], "Element not in Formula upon decrement");
196 elementCounts[Z-1]--; // atomic numbers start at 1
197 // might need to update number of elements
198 if(!elementCounts[Z-1]){
199 numElements--;
200 // resize the Array if this was at the last position
201 if(Z==elementCounts.size()){
202 // find the first element from the back that is not equal to zero
203 set_t::reverse_iterator riter = find_if(elementCounts.rbegin(),
204 elementCounts.rend(),
205 bind1st(not_equal_to<mapped_type>(),0));
206 // see how many elements are in this range
207 set_t::reverse_iterator::difference_type diff = riter - elementCounts.rbegin();
208 elementCounts.resize(elementCounts.size()-diff);
209 }
210 }
211}
212
213void Formula::operator-=(const string &shorthand){
214 const element * element = World::getInstance().getPeriode()->FindElement(shorthand);
215 operator-=(element);
216}
217
218void Formula::addElements(const element *element,unsigned int count){
219 ASSERT(element,"Invalid pointer in Formula::addElements(element*)");
220 addElements(element->getAtomicNumber(),count);
221}
222
223void Formula::addElements(atomicNumber_t Z,unsigned int count){
224 if(count==0) return;
225 ASSERT(Z>0,"Invalid atomic Number");
226 ASSERT(World::getInstance().getPeriode()->FindElement(Z),"No Element with this number in Periodentafel");
227 elementCounts.resize(max<atomicNumber_t>(Z,elementCounts.size()),0); // No-op when we already have the right size
228 // might need to update number of elements
229 if(!elementCounts[Z-1]){
230 numElements++;
231 }
232 elementCounts[Z-1]+=count;
233}
234
235void Formula::addElements(const string &shorthand,unsigned int count){
236 const element * element = World::getInstance().getPeriode()->FindElement(shorthand);
237 addElements(element,count);
238}
239
240void Formula::addFormula(const Formula &formula,unsigned int n){
241 for(Formula::const_iterator iter=formula.begin();iter!=formula.end();++iter){
242 this->addElements(iter->first,iter->second*n);
243 }
244}
245
246enumeration<Formula::key_type> Formula::enumerateElements() const{
247 enumeration<key_type> res(1);
248 for(Formula::const_iterator iter=begin();iter!=end();++iter){
249 res.add(iter->first);
250 }
251 return res;
252}
253
254const unsigned int Formula::operator[](const element *element) const{
255 ASSERT(element,"Invalid pointer in access of Formula");
256 return operator[](element->getAtomicNumber());
257}
258
259const unsigned int Formula::operator[](atomicNumber_t Z) const{
260 ASSERT(Z>0,"Invalid atomic Number");
261 ASSERT(World::getInstance().getPeriode()->FindElement(Z),"No Element with this number in Periodentafel");
262 if(elementCounts.size()<Z)
263 return 0;
264 return elementCounts[Z-1]; // atomic numbers start at 1
265}
266
267const unsigned int Formula::operator[](std::string shorthand) const{
268 const element * element = World::getInstance().getPeriode()->FindElement(shorthand);
269 return operator[](element);
270}
271
272bool Formula::operator==(const Formula &rhs) const{
273 // quick check... number of elements used
274 if(numElements != rhs.numElements){
275 return false;
276 }
277 // second quick check, size of vectors (== last element in formula)
278 if(elementCounts.size()!=rhs.elementCounts.size()){
279 return false;
280 }
281 // slow check: all elements
282 // direct access to internal structure means all element-counts have to be compared
283 // this avoids access to periodentafel to find elements though and is probably faster
284 // in total
285 return equal(elementCounts.begin(),
286 elementCounts.end(),
287 rhs.elementCounts.begin());
288}
289
290bool Formula::operator!=(const Formula &rhs) const{
291 return !operator==(rhs);
292}
293
294Formula::iterator Formula::begin(){
295 return iterator(elementCounts,0);
296}
297Formula::const_iterator Formula::begin() const{
298 // this is the only place where this is needed, so this is better than making it mutable
299 return const_iterator(const_cast<set_t&>(elementCounts),0);
300}
301Formula::iterator Formula::end(){
302 return iterator(elementCounts);
303}
304Formula::const_iterator Formula::end() const{
305 // this is the only place where this is needed, so this is better than making it mutable
306 return const_iterator(const_cast<set_t&>(elementCounts));
307}
308
309void Formula::clear(){
310 elementCounts.clear();
311 numElements = 0;
312}
313
314/**************** Iterator structure ********************/
315
316template <class result_type>
317Formula::_iterator<result_type>::_iterator(set_t &_set) :
318 set(&_set)
319{
320 pos=set->size();
321}
322
323template <class result_type>
324Formula::_iterator<result_type>::_iterator(set_t &_set,size_t _pos) :
325 set(&_set),pos(_pos)
326{
327 ASSERT(pos<=set->size(),"invalid position in iterator construction");
328 while(pos<set->size() && (*set)[pos]==0) ++pos;
329}
330
331template <class result_type>
332Formula::_iterator<result_type>::_iterator(const _iterator &rhs) :
333 set(rhs.set),pos(rhs.pos)
334{}
335
336template <class result_type>
337Formula::_iterator<result_type>::~_iterator(){}
338
339template <class result_type>
340Formula::_iterator<result_type>&
341Formula::_iterator<result_type>::operator=(const _iterator<result_type> &rhs){
342 set=rhs.set;
343 pos=rhs.pos;
344 return *this;
345}
346
347template <class result_type>
348bool
349Formula::_iterator<result_type>::operator==(const _iterator<result_type> &rhs){
350 return set==rhs.set && pos==rhs.pos;
351}
352
353template <class result_type>
354bool
355Formula::_iterator<result_type>::operator!=(const _iterator<result_type> &rhs){
356 return !operator==(rhs);
357}
358
359template <class result_type>
360Formula::_iterator<result_type>
361Formula::_iterator<result_type>::operator++(){
362 ASSERT(pos!=set->size(),"Incrementing Formula::iterator beyond end");
363 pos++;
364 while(pos<set->size() && (*set)[pos]==0) ++pos;
365 return *this;
366}
367
368template <class result_type>
369Formula::_iterator<result_type>
370Formula::_iterator<result_type>::operator++(int){
371 Formula::_iterator<result_type> retval = *this;
372 ++(*this);
373 return retval;
374}
375
376template <class result_type>
377Formula::_iterator<result_type>
378Formula::_iterator<result_type>::operator--(){
379 ASSERT(pos!=0,"Decrementing Formula::iterator beyond begin");
380 pos--;
381 while(pos>0 && (*set)[pos]==0) --pos;
382 return *this;
383}
384
385template <class result_type>
386Formula::_iterator<result_type>
387Formula::_iterator<result_type>::operator--(int){
388 Formula::_iterator<result_type> retval = *this;
389 --(*this);
390 return retval;
391}
392
393template <class result_type>
394result_type
395Formula::_iterator<result_type>::operator*(){
396 const element *element = World::getInstance().getPeriode()->FindElement(pos+1);
397 ASSERT(element,"Element with position of iterator not found");
398 return make_pair(element,(*set)[pos]);
399}
400
401template <class result_type>
402result_type*
403Formula::_iterator<result_type>::operator->(){
404 // no one can keep this value around, so a static is ok to avoid temporaries
405 static value_type value=make_pair(reinterpret_cast<element*>(0),0); // no default constructor for std::pair
406 const element *element = World::getInstance().getPeriode()->FindElement(pos+1);
407 ASSERT(element,"Element with position of iterator not found");
408 value = make_pair(element,(*set)[pos]);
409 return &value;
410}
411
412// explicit instantiation of all iterator template methods
413// this is quite ugly, but there is no better way unless we expose iterator implementation
414
415// instantiate Formula::iterator
416template Formula::iterator::_iterator(set_t&);
417template Formula::iterator::_iterator(set_t&,size_t);
418template Formula::iterator::_iterator(const Formula::iterator&);
419template Formula::iterator::~_iterator();
420template Formula::iterator &Formula::iterator::operator=(const Formula::iterator&);
421template bool Formula::iterator::operator==(const Formula::iterator&);
422template bool Formula::iterator::operator!=(const Formula::iterator&);
423template Formula::iterator Formula::iterator::operator++();
424template Formula::iterator Formula::iterator::operator++(int);
425template Formula::iterator Formula::iterator::operator--();
426template Formula::iterator Formula::iterator::operator--(int);
427template Formula::value_type Formula::iterator::operator*();
428template Formula::value_type *Formula::iterator::operator->();
429
430// instantiate Formula::const_iterator
431template Formula::const_iterator::_iterator(set_t&);
432template Formula::const_iterator::_iterator(set_t&,size_t);
433template Formula::const_iterator::_iterator(const Formula::const_iterator&);
434template Formula::const_iterator::~_iterator();
435template Formula::const_iterator &Formula::const_iterator::operator=(const Formula::const_iterator&);
436template bool Formula::const_iterator::operator==(const Formula::const_iterator&);
437template bool Formula::const_iterator::operator!=(const Formula::const_iterator&);
438template Formula::const_iterator Formula::const_iterator::operator++();
439template Formula::Formula::const_iterator Formula::const_iterator::operator++(int);
440template Formula::Formula::const_iterator Formula::const_iterator::operator--();
441template Formula::Formula::const_iterator Formula::const_iterator::operator--(int);
442template const Formula::value_type Formula::const_iterator::operator*();
443template const Formula::value_type *Formula::const_iterator::operator->();
444
445/********************** I/O of Formulas ************************************************/
446
447std::ostream &operator<<(std::ostream &ost,const Formula &formula){
448 ost << formula.toString();
449 return ost;
450}
Note: See TracBrowser for help on using the repository browser.