source: src/molecule.cpp@ 4464ef

Action_Thermostats Add_AtomRandomPerturbation Add_RotateAroundBondAction Add_SelectAtomByNameAction Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_StructOpt_integration_tests Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_ChronosMutex Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion GeometryObjects Gui_displays_atomic_force_velocity IndependentFragmentGrids_IntegrationTest JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks RotateToPrincipalAxisSystem_UndoRedo StoppableMakroAction Subpackage_levmar Subpackage_vmg ThirdParty_MPQC_rebuilt_buildsystem TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps Ubuntu_1604_changes stable
Last change on this file since 4464ef was 1c9eb8, checked in by Frederik Heber <heber@…>, 9 years ago

FIX: molecule::SetNameFromFilename did not notify about MoleculeNameChanged.

  • Property mode set to 100755
File size: 46.8 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 * Copyright (C) 2013 Frederik Heber. All rights reserved.
6 *
7 *
8 * This file is part of MoleCuilder.
9 *
10 * MoleCuilder is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * MoleCuilder is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
22 */
23
24/** \file molecules.cpp
25 *
26 * Functions for the class molecule.
27 *
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 <algorithm>
38#include <boost/assign.hpp>
39#include <boost/bind.hpp>
40#include <boost/foreach.hpp>
41#include <cstring>
42
43#include <gsl/gsl_inline.h>
44#include <gsl/gsl_heapsort.h>
45
46#include "molecule.hpp"
47
48#include "Atom/atom.hpp"
49#include "Bond/bond.hpp"
50#include "Box.hpp"
51#include "CodePatterns/enumeration.hpp"
52#include "CodePatterns/Log.hpp"
53#include "CodePatterns/Observer/Observable.hpp"
54#include "CodePatterns/Observer/Notification.hpp"
55#include "config.hpp"
56#include "Descriptors/AtomIdDescriptor.hpp"
57#include "Element/element.hpp"
58#include "Graph/BondGraph.hpp"
59#include "LinearAlgebra/Exceptions.hpp"
60#include "LinearAlgebra/leastsquaremin.hpp"
61#include "LinearAlgebra/Plane.hpp"
62#include "LinearAlgebra/RealSpaceMatrix.hpp"
63#include "LinearAlgebra/Vector.hpp"
64#include "LinkedCell/linkedcell.hpp"
65#include "MoleculeObserver.hpp"
66#include "IdPool_impl.hpp"
67#include "Shapes/BaseShapes.hpp"
68#include "Tesselation/tesselation.hpp"
69#include "World.hpp"
70#include "WorldTime.hpp"
71
72using namespace boost::assign;
73
74// static entities
75static Observable::channels_t getAtomPositionsChannels()
76{
77 Observable::channels_t channels;
78 channels += molecule::AtomInserted, molecule::AtomRemoved, molecule::AtomMoved;
79 return channels;
80}
81
82/************************************* Functions for class molecule *********************************/
83
84/** Constructor of class molecule.
85 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
86 */
87molecule::molecule() :
88 Observable("molecule"),
89 MDSteps(0),
90 NoNonBonds(0),
91 NoCyclicBonds(0),
92 ActiveFlag(false),
93 IndexNr(-1),
94 NoNonHydrogen(0),
95 BondCount(0),
96 atomIdPool(1, 20, 100),
97 BoundingBoxSweepingAxis(std::vector<AtomDistanceMap_t>(NDIM)),
98 _lastchangedatomid(-1),
99 last_atom(0),
100 molcenter(zeroVec),
101 selected(false)
102{
103 // add specific channels
104 Channels *OurChannel = new Channels;
105 Observable::insertNotificationChannel( std::make_pair( static_cast<Observable *>(this), OurChannel) );
106 for (size_t type = 0; type < (size_t)NotificationType_MAX; ++type)
107 OurChannel->addChannel(type);
108
109 // cannot initialize in initializer body as then channels have not been setup yet
110 BoundingBox.reset(
111 new Cacheable<BoundingBoxInfo>(
112 this, boost::bind(&molecule::updateBoundingBox, this), "molecule_BoundingBox", getAtomPositionsChannels()));
113 MoleculeCenter.reset(
114 new Cacheable<Vector>(
115 this, boost::bind(&molecule::updateMoleculeCenter, this), "molecule_center", getAtomPositionsChannels()));
116
117 strcpy(name,World::getInstance().getDefaultName().c_str());
118
119 // inform MoleculeObserver about new molecule
120 MoleculeObserver::getInstance().Inserted(this);
121}
122
123molecule *NewMolecule(){
124 return new molecule();
125}
126
127/** Destructor of class molecule.
128 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
129 */
130molecule::~molecule()
131{
132 CleanupMolecule();
133
134 // inform MoleculeObserver about removed molecule
135 MoleculeObserver::getInstance().Removed(this);
136};
137
138
139void DeleteMolecule(molecule *mol){
140 delete mol;
141}
142
143// getter and setter
144const std::string molecule::getName() const{
145 return std::string(name);
146}
147
148int molecule::getAtomCount() const{
149 return atomIds.size();
150}
151
152void molecule::setName(const std::string _name){
153 OBSERVE;
154 NOTIFY(MoleculeNameChanged);
155 cout << "Set name of molecule " << getId() << " to " << _name << endl;
156 strncpy(name,_name.c_str(),MAXSTRINGSIZE);
157}
158
159void molecule::InsertLocalToGlobalId(atom * const pointer)
160{
161#ifndef NDEBUG
162 std::pair< LocalToGlobalId_t::iterator, bool > inserter =
163#endif
164 LocalToGlobalId.insert( std::make_pair(pointer->getNr(), pointer) );
165 ASSERT( inserter.second,
166 "molecule::AddAtom() - local number "+toString(pointer->getNr())+" appears twice.");
167}
168
169bool molecule::changeAtomNr(int oldNr, int newNr, atom* target){
170 OBSERVE;
171 if(atomIdPool.reserveId(newNr)){
172 NOTIFY(AtomNrChanged);
173 if (oldNr != -1) // -1 is reserved and indicates no number
174 atomIdPool.releaseId(oldNr);
175 LocalToGlobalId.erase(oldNr);
176 ASSERT (target,
177 "molecule::changeAtomNr() - given target is NULL, cannot set Nr or name.");
178 target->setNr(newNr);
179 _lastchangedatomid = target->getId();
180 InsertLocalToGlobalId(target);
181 setAtomName(target);
182 return true;
183 } else{
184 return false;
185 }
186}
187
188bool molecule::changeAtomId(int oldId, int newId)
189{
190 OBSERVE;
191 if ((!atomIds.contains( oldId )) || (atomIds.contains( newId )))
192 return false;
193 atomIds.erase( oldId );
194 atomIds.insert( newId );
195 // also update BoundingBoxSweepingAxis
196 for (int i=0;i<NDIM;++i) {
197 AtomDistanceMap_t::left_iterator iter = BoundingBoxSweepingAxis[i].left.find(oldId);
198 ASSERT(iter != BoundingBoxSweepingAxis[i].left.end(),
199 "molecule::changeAtomId() - could not find atom "+toString(oldId)
200 +" in BoundingBoxSweepingAxis.");
201 const double component = iter->second;
202 BoundingBoxSweepingAxis[i].left.erase(iter);
203 BoundingBoxSweepingAxis[i].left.insert( std::make_pair(newId, component) );
204 }
205 return true;
206}
207
208bool molecule::changeId(moleculeId_t newId){
209 // first we move ourselves in the world
210 // the world lets us know if that succeeded
211 if(World::getInstance().changeMoleculeId(id,newId,this)){
212 OBSERVE;
213 NOTIFY(IndexChanged);
214 id = newId;
215 return true;
216 }
217 else{
218 return false;
219 }
220}
221
222
223moleculeId_t molecule::getId() const {
224 return id;
225}
226
227void molecule::setId(moleculeId_t _id){
228 id =_id;
229}
230
231const Formula &molecule::getFormula() const {
232 return formula;
233}
234
235unsigned int molecule::getElementCount() const{
236 return formula.getElementCount();
237}
238
239bool molecule::hasElement(const element *element) const{
240 return formula.hasElement(element);
241}
242
243bool molecule::hasElement(atomicNumber_t Z) const{
244 return formula.hasElement(Z);
245}
246
247bool molecule::hasElement(const string &shorthand) const{
248 return formula.hasElement(shorthand);
249}
250
251/************************** Access to the List of Atoms ****************/
252
253molecule::const_iterator molecule::erase( const_iterator loc )
254{
255 OBSERVE;
256 const_iterator iter = loc;
257 ++iter;
258 atom * const _atom = const_cast<atom *>(*loc);
259 {
260 _lastchangedatomid = _atom->getId();
261 NOTIFY(AtomRemoved);
262 }
263 atomIds.erase( _atom->getId() );
264 {
265 BoundingBoxInfo oldinfo = updateBoundingBox();
266 for (int i=0;i<NDIM;++i)
267 BoundingBoxSweepingAxis[i].left.erase( _atom->getId() );
268 BoundingBoxInfo newinfo = updateBoundingBox();
269 if (oldinfo != newinfo)
270 NOTIFY(BoundingBoxChanged);
271 }
272 {
273 molcenter -= _atom->getPosition();
274 }
275 {
276 NOTIFY(AtomNrChanged);
277 atomIdPool.releaseId(_atom->getNr());
278 LocalToGlobalId.erase(_atom->getNr());
279 _atom->setNr(-1);
280 }
281 _atom->removeFromMolecule();
282 return iter;
283}
284
285molecule::const_iterator molecule::erase( atom * key )
286{
287 const_iterator iter = const_cast<const molecule &>(*this).find(key);
288 if (iter != const_cast<const molecule &>(*this).end())
289 return erase(iter);
290 else
291 return iter;
292}
293
294pair<molecule::iterator,bool> molecule::insert ( atom * const key )
295{
296 OBSERVE;
297 NOTIFY(AtomInserted);
298 _lastchangedatomid = key->getId();
299 std::pair<iterator,bool> res = atomIds.insert(key->getId());
300 if (res.second) { // push atom if went well
301 {
302 BoundingBoxInfo oldinfo = updateBoundingBox();
303 for (int i=0;i<NDIM;++i)
304 BoundingBoxSweepingAxis[i].left.insert( std::make_pair(key->getId(), key->getPosition()[i]));
305 BoundingBoxInfo newinfo = updateBoundingBox();
306 if (oldinfo != newinfo)
307 NOTIFY(BoundingBoxChanged);
308 }
309 {
310 molcenter += key->getPosition();
311 }
312 NOTIFY(AtomNrChanged);
313 key->setNr(atomIdPool.getNextId());
314 InsertLocalToGlobalId(key);
315 setAtomName(key);
316 return res;
317 } else {
318 return pair<iterator,bool>(end(),res.second);
319 }
320}
321
322void molecule::setAtomName(atom *_atom) const
323{
324 std::stringstream sstr;
325 sstr << _atom->getType()->getSymbol() << _atom->getNr();
326 _atom->setName(sstr.str());
327}
328
329World::AtomComposite molecule::getAtomSet()
330{
331 World::AtomComposite vector_of_atoms;
332 for (molecule::iterator iter = begin(); iter != end(); ++iter)
333 vector_of_atoms.push_back(*iter);
334 return vector_of_atoms;
335}
336
337World::ConstAtomComposite molecule::getAtomSet() const
338{
339 World::ConstAtomComposite vector_of_atoms;
340 for (molecule::const_iterator iter = begin(); iter != end(); ++iter)
341 vector_of_atoms.push_back(*iter);
342 return vector_of_atoms;
343}
344
345/** Adds given atom \a *pointer from molecule list.
346 * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount
347 * \param *pointer allocated and set atom
348 * \return true - succeeded, false - atom not found in list
349 */
350bool molecule::AddAtom(atom *pointer)
351{
352 if (pointer != NULL) {
353 // molecule::insert() is called by setMolecule()
354 pointer->setMolecule(this);
355 }
356 return true;
357};
358
359/** Adds a copy of the given atom \a *pointer from molecule list.
360 * Increases molecule::last_atom and gives last number to added atom.
361 * \param *pointer allocated and set atom
362 * \return pointer to the newly added atom
363 */
364atom * molecule::AddCopyAtom(atom *pointer)
365{
366 atom *retval = NULL;
367 if (pointer != NULL) {
368 atom *walker = pointer->clone();
369 AddAtom(walker);
370 retval=walker;
371 }
372 return retval;
373};
374
375/** Adds a Hydrogen atom in replacement for the given atom \a *partner in bond with a *origin.
376 * Here, we have to distinguish between single, double or triple bonds as stated by \a BondDegree, that each demand
377 * a different scheme when adding \a *replacement atom for the given one.
378 * -# Single Bond: Simply add new atom with bond distance rescaled to typical hydrogen one
379 * -# Double Bond: Here, we need the **BondList of the \a *origin atom, by scanning for the other bonds instead of
380 * *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalvector().
381 * The orthonormal vector to this plane along with the vector in *Bond direction determines the plane the two
382 * replacing hydrogens shall lie in. Now, all remains to do is take the usual hydrogen double bond angle for the
383 * element of *origin and form the sin/cos admixture of both plane vectors for the new coordinates of the two
384 * hydrogens forming this angle with *origin.
385 * -# Triple Bond: The idea is to set up a tetraoid (C1-H1-H2-H3) (however the lengths \f$b\f$ of the sides of the base
386 * triangle formed by the to be added hydrogens are not equal to the typical bond distance \f$l\f$ but have to be
387 * determined from the typical angle \f$\alpha\f$ for a hydrogen triple connected to the element of *origin):
388 * We have the height \f$d\f$ as the vector in *Bond direction (from triangle C1-H1-H2).
389 * \f[ h = l \cdot \cos{\left (\frac{\alpha}{2} \right )} \qquad b = 2l \cdot \sin{\left (\frac{\alpha}{2} \right)} \quad \rightarrow \quad d = l \cdot \sqrt{\cos^2{\left (\frac{\alpha}{2} \right)}-\frac{1}{3}\cdot\sin^2{\left (\frac{\alpha}{2}\right )}}
390 * \f]
391 * vector::GetNormalvector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalvector creates
392 * the third one from the former two vectors. The latter ones form the plane of the base triangle mentioned above.
393 * The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
394 * the median lines in an isosceles triangle meet in the center point with a ratio 2:1.
395 * \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2}
396 * \f]
397 * as the coordination of all three atoms in the coordinate system of these three vectors:
398 * \f$\pmatrix{d & f & 0}\f$, \f$\pmatrix{d & -0.5 \cdot f & g}\f$ and \f$\pmatrix{d & -0.5 \cdot f & -g}\f$.
399 *
400 * \param *out output stream for debugging
401 * \param *Bond pointer to bond between \a *origin and \a *replacement
402 * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin)
403 * \param *origin pointer to atom which acts as the origin for scaling the added hydrogen to correct bond length
404 * \param *replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
405 * \param isAngstroem whether the coordination of the given atoms is in AtomicLength (false) or Angstrom(true)
406 * \return number of atoms added, if < bond::BondDegree then something went wrong
407 * \todo double and triple bonds splitting (always use the tetraeder angle!)
408 */
409//bool molecule::AddHydrogenReplacementAtom(bond::ptr TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
410//{
411//// Info info(__func__);
412// bool AllWentWell = true; // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
413// double bondlength; // bond length of the bond to be replaced/cut
414// double bondangle; // bond angle of the bond to be replaced/cut
415// double BondRescale; // rescale value for the hydrogen bond length
416// bond::ptr FirstBond;
417// bond::ptr SecondBond; // Other bonds in double bond case to determine "other" plane
418// atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
419// double b,l,d,f,g, alpha, factors[NDIM]; // hold temporary values in triple bond case for coordination determination
420// Vector Orthovector1, Orthovector2; // temporary vectors in coordination construction
421// Vector InBondvector; // vector in direction of *Bond
422// const RealSpaceMatrix &matrix = World::getInstance().getDomain().getM();
423// bond::ptr Binder;
424//
425// // create vector in direction of bond
426// InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
427// bondlength = InBondvector.Norm();
428//
429// // is greater than typical bond distance? Then we have to correct periodically
430// // the problem is not the H being out of the box, but InBondvector have the wrong direction
431// // due to TopReplacement or Origin being on the wrong side!
432// const BondGraph * const BG = World::getInstance().getBondGraph();
433// const range<double> MinMaxBondDistance(
434// BG->getMinMaxDistance(TopOrigin,TopReplacement));
435// if (!MinMaxBondDistance.isInRange(bondlength)) {
436//// LOG(4, "InBondvector is: " << InBondvector << ".");
437// Orthovector1.Zero();
438// for (int i=NDIM;i--;) {
439// l = TopReplacement->at(i) - TopOrigin->at(i);
440// if (fabs(l) > MinMaxBondDistance.last) { // is component greater than bond distance (check against min not useful here)
441// Orthovector1[i] = (l < 0) ? -1. : +1.;
442// } // (signs are correct, was tested!)
443// }
444// Orthovector1 *= matrix;
445// InBondvector -= Orthovector1; // subtract just the additional translation
446// bondlength = InBondvector.Norm();
447//// LOG(4, "INFO: Corrected InBondvector is now: " << InBondvector << ".");
448// } // periodic correction finished
449//
450// InBondvector.Normalize();
451// // get typical bond length and store as scale factor for later
452// ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
453// BondRescale = TopOrigin->getType()->getHBondDistance(TopBond->getDegree()-1);
454// if (BondRescale == -1) {
455// ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
456// return false;
457// BondRescale = bondlength;
458// } else {
459// if (!IsAngstroem)
460// BondRescale /= (1.*AtomicLengthToAngstroem);
461// }
462//
463// // discern single, double and triple bonds
464// switch(TopBond->getDegree()) {
465// case 1:
466// FirstOtherAtom = World::getInstance().createAtom(); // new atom
467// FirstOtherAtom->setType(1); // element is Hydrogen
468// FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
469// FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
470// if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
471// FirstOtherAtom->father = TopReplacement;
472// BondRescale = bondlength;
473// } else {
474// FirstOtherAtom->father = NULL; // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
475// }
476// InBondvector *= BondRescale; // rescale the distance vector to Hydrogen bond length
477// FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
478// AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
479//// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
480// Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
481// Binder->Cyclic = false;
482// Binder->Type = GraphEdge::TreeEdge;
483// break;
484// case 2:
485// {
486// // determine two other bonds (warning if there are more than two other) plus valence sanity check
487// const BondList& ListOfBonds = TopOrigin->getListOfBonds();
488// for (BondList::const_iterator Runner = ListOfBonds.begin();
489// Runner != ListOfBonds.end();
490// ++Runner) {
491// if ((*Runner) != TopBond) {
492// if (FirstBond == NULL) {
493// FirstBond = (*Runner);
494// FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
495// } else if (SecondBond == NULL) {
496// SecondBond = (*Runner);
497// SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
498// } else {
499// ELOG(2, "Detected more than four bonds for atom " << TopOrigin->getName());
500// }
501// }
502// }
503// }
504// if (SecondOtherAtom == NULL) { // then we have an atom with valence four, but only 3 bonds: one to replace and one which is TopBond (third is FirstBond)
505// SecondBond = TopBond;
506// SecondOtherAtom = TopReplacement;
507// }
508// if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
509//// LOG(3, "Regarding the double bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << TopOrigin->Name << " to determine orthogonal plane.");
510//
511// // determine the plane of these two with the *origin
512// try {
513// Orthovector1 = Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
514// }
515// catch(LinearDependenceException &excp){
516// LOG(0, boost::diagnostic_information(excp));
517// // TODO: figure out what to do with the Orthovector in this case
518// AllWentWell = false;
519// }
520// } else {
521// Orthovector1.GetOneNormalVector(InBondvector);
522// }
523// //LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
524// // orthogonal vector and bond vector between origin and replacement form the new plane
525// Orthovector1.MakeNormalTo(InBondvector);
526// Orthovector1.Normalize();
527// //LOG(3, "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << ".");
528//
529// // create the two Hydrogens ...
530// FirstOtherAtom = World::getInstance().createAtom();
531// SecondOtherAtom = World::getInstance().createAtom();
532// FirstOtherAtom->setType(1);
533// SecondOtherAtom->setType(1);
534// FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
535// FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
536// SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
537// SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
538// FirstOtherAtom->father = NULL; // we are just an added hydrogen with no father
539// SecondOtherAtom->father = NULL; // we are just an added hydrogen with no father
540// bondangle = TopOrigin->getType()->getHBondAngle(1);
541// if (bondangle == -1) {
542// ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
543// return false;
544// bondangle = 0;
545// }
546// bondangle *= M_PI/180./2.;
547//// LOG(3, "INFO: ReScaleCheck: InBondvector " << InBondvector << ", " << Orthovector1 << ".");
548//// LOG(3, "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle));
549// FirstOtherAtom->Zero();
550// SecondOtherAtom->Zero();
551// for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
552// FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
553// SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
554// }
555// FirstOtherAtom->Scale(BondRescale); // rescale by correct BondDistance
556// SecondOtherAtom->Scale(BondRescale);
557// //LOG(3, "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << ".");
558// *FirstOtherAtom += TopOrigin->getPosition();
559// *SecondOtherAtom += TopOrigin->getPosition();
560// // ... and add to molecule
561// AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
562// AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
563//// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
564//// LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
565// Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
566// Binder->Cyclic = false;
567// Binder->Type = GraphEdge::TreeEdge;
568// Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
569// Binder->Cyclic = false;
570// Binder->Type = GraphEdge::TreeEdge;
571// break;
572// case 3:
573// // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
574// FirstOtherAtom = World::getInstance().createAtom();
575// SecondOtherAtom = World::getInstance().createAtom();
576// ThirdOtherAtom = World::getInstance().createAtom();
577// FirstOtherAtom->setType(1);
578// SecondOtherAtom->setType(1);
579// ThirdOtherAtom->setType(1);
580// FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
581// FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
582// SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
583// SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
584// ThirdOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
585// ThirdOtherAtom->setFixedIon(TopReplacement->getFixedIon());
586// FirstOtherAtom->father = NULL; // we are just an added hydrogen with no father
587// SecondOtherAtom->father = NULL; // we are just an added hydrogen with no father
588// ThirdOtherAtom->father = NULL; // we are just an added hydrogen with no father
589//
590// // we need to vectors orthonormal the InBondvector
591// AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
592//// LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
593// try{
594// Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
595// }
596// catch(LinearDependenceException &excp) {
597// LOG(0, boost::diagnostic_information(excp));
598// AllWentWell = false;
599// }
600//// LOG(3, "INFO: Orthovector2: " << Orthovector2 << ".")
601//
602// // create correct coordination for the three atoms
603// alpha = (TopOrigin->getType()->getHBondAngle(2))/180.*M_PI/2.; // retrieve triple bond angle from database
604// l = BondRescale; // desired bond length
605// b = 2.*l*sin(alpha); // base length of isosceles triangle
606// d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.); // length for InBondvector
607// f = b/sqrt(3.); // length for Orthvector1
608// g = b/2.; // length for Orthvector2
609//// LOG(3, "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", ");
610//// LOG(3, "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g));
611// factors[0] = d;
612// factors[1] = f;
613// factors[2] = 0.;
614// FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
615// factors[1] = -0.5*f;
616// factors[2] = g;
617// SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
618// factors[2] = -g;
619// ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
620//
621// // rescale each to correct BondDistance
622//// FirstOtherAtom->x.Scale(&BondRescale);
623//// SecondOtherAtom->x.Scale(&BondRescale);
624//// ThirdOtherAtom->x.Scale(&BondRescale);
625//
626// // and relative to *origin atom
627// *FirstOtherAtom += TopOrigin->getPosition();
628// *SecondOtherAtom += TopOrigin->getPosition();
629// *ThirdOtherAtom += TopOrigin->getPosition();
630//
631// // ... and add to molecule
632// AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
633// AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
634// AllWentWell = AllWentWell && AddAtom(ThirdOtherAtom);
635//// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
636//// LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
637//// LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
638// Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
639// Binder->Cyclic = false;
640// Binder->Type = GraphEdge::TreeEdge;
641// Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
642// Binder->Cyclic = false;
643// Binder->Type = GraphEdge::TreeEdge;
644// Binder = AddBond(BottomOrigin, ThirdOtherAtom, 1);
645// Binder->Cyclic = false;
646// Binder->Type = GraphEdge::TreeEdge;
647// break;
648// default:
649// ELOG(1, "BondDegree does not state single, double or triple bond!");
650// AllWentWell = false;
651// break;
652// }
653//
654// return AllWentWell;
655//};
656
657/** Creates a copy of this molecule.
658 * \param offset translation Vector for the new molecule relative to old one
659 * \return copy of molecule
660 */
661molecule *molecule::CopyMolecule(const Vector &offset)
662{
663 molecule *copy = World::getInstance().createMolecule();
664
665 // copy all atoms
666 std::map< const atom *, atom *> FatherFinder;
667 for (iterator iter = begin(); iter != end(); ++iter) {
668 atom * const copy_atom = copy->AddCopyAtom(*iter);
669 copy_atom->setPosition(copy_atom->getPosition() + offset);
670 FatherFinder.insert( std::make_pair( *iter, copy_atom ) );
671 }
672
673 // copy all bonds
674 for(const_iterator AtomRunner = const_cast<const molecule &>(*this).begin();
675 AtomRunner != const_cast<const molecule &>(*this).end();
676 ++AtomRunner) {
677 const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
678 for(BondList::const_iterator BondRunner = ListOfBonds.begin();
679 BondRunner != ListOfBonds.end();
680 ++BondRunner)
681 if ((*BondRunner)->leftatom == *AtomRunner) {
682 bond::ptr Binder = (*BondRunner);
683 // get the pendant atoms of current bond in the copy molecule
684 ASSERT(FatherFinder.count(Binder->leftatom),
685 "molecule::CopyMolecule() - No copy of original left atom "
686 +toString(Binder->leftatom)+" for bond copy found");
687 ASSERT(FatherFinder.count(Binder->rightatom),
688 "molecule::CopyMolecule() - No copy of original right atom "
689 +toString(Binder->rightatom)+" for bond copy found");
690 atom * const LeftAtom = FatherFinder[Binder->leftatom];
691 atom * const RightAtom = FatherFinder[Binder->rightatom];
692
693 bond::ptr const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->getDegree());
694 NewBond->Cyclic = Binder->Cyclic;
695 if (Binder->Cyclic)
696 copy->NoCyclicBonds++;
697 NewBond->Type = Binder->Type;
698 }
699 }
700 // correct fathers
701 //for_each(begin(),end(),mem_fun(&atom::CorrectFather));
702
703 return copy;
704};
705
706
707/** Destroys all atoms inside this molecule.
708 */
709void removeAtomsinMolecule(molecule *&_mol)
710{
711 // copy list of atoms from molecule as it will be changed
712 std::vector<atom *> atoms;
713 atoms.resize(_mol->getAtomCount(), NULL);
714 std::copy(_mol->begin(), _mol->end(), atoms.begin());
715 // remove each atom from world
716 for(std::vector<atom *>::iterator AtomRunner = atoms.begin();
717 AtomRunner != atoms.end(); ++AtomRunner)
718 World::getInstance().destroyAtom(*AtomRunner);
719 // make sure that pointer os not usable
720 _mol = NULL;
721};
722
723
724/**
725 * Copies all atoms of a molecule which are within the defined parallelepiped.
726 *
727 * @param offest for the origin of the parallelepiped
728 * @param three vectors forming the matrix that defines the shape of the parallelpiped
729 */
730molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) {
731 molecule *copy = World::getInstance().createMolecule();
732
733 // copy all atoms
734 std::map< const atom *, atom *> FatherFinder;
735 for (iterator iter = begin(); iter != end(); ++iter) {
736 if (region.isInside((*iter)->getPosition())) {
737 atom * const copy_atom = copy->AddCopyAtom(*iter);
738 FatherFinder.insert( std::make_pair( *iter, copy_atom ) );
739 }
740 }
741
742 // copy all bonds
743 for(molecule::const_iterator AtomRunner = const_cast<const molecule &>(*this).begin();
744 AtomRunner != const_cast<const molecule &>(*this).end();
745 ++AtomRunner) {
746 const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
747 for(BondList::const_iterator BondRunner = ListOfBonds.begin();
748 BondRunner != ListOfBonds.end();
749 ++BondRunner)
750 if ((*BondRunner)->leftatom == *AtomRunner) {
751 bond::ptr Binder = (*BondRunner);
752 if ((FatherFinder.count(Binder->leftatom))
753 && (FatherFinder.count(Binder->rightatom))) {
754 // if copy present, then it must be from subregion
755 atom * const LeftAtom = FatherFinder[Binder->leftatom];
756 atom * const RightAtom = FatherFinder[Binder->rightatom];
757
758 bond::ptr const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->getDegree());
759 NewBond->Cyclic = Binder->Cyclic;
760 if (Binder->Cyclic)
761 copy->NoCyclicBonds++;
762 NewBond->Type = Binder->Type;
763 }
764 }
765 }
766 // correct fathers
767 //for_each(begin(),end(),mem_fun(&atom::CorrectFather));
768
769 //TODO: copy->BuildInducedSubgraph(this);
770
771 return copy;
772}
773
774/** Adds a bond to a the molecule specified by two atoms, \a *first and \a *second.
775 * Also updates molecule::BondCount and molecule::NoNonBonds.
776 * \param *first first atom in bond
777 * \param *second atom in bond
778 * \return pointer to bond or NULL on failure
779 */
780bond::ptr molecule::AddBond(atom *atom1, atom *atom2, int degree)
781{
782 bond::ptr Binder;
783
784 // some checks to make sure we are able to create the bond
785 ASSERT(atom1,
786 "molecule::AddBond() - First atom "+toString(atom1)
787 +" is not a invalid pointer");
788 ASSERT(atom2,
789 "molecule::AddBond() - Second atom "+toString(atom2)
790 +" is not a invalid pointer");
791 ASSERT(isInMolecule(atom1),
792 "molecule::AddBond() - First atom "+toString(atom1)
793 +" is not part of molecule");
794 ASSERT(isInMolecule(atom2),
795 "molecule::AddBond() - Second atom "+toString(atom2)
796 +" is not part of molecule");
797
798 Binder.reset(new bond(atom1, atom2, degree));
799 atom1->RegisterBond(WorldTime::getTime(), Binder);
800 atom2->RegisterBond(WorldTime::getTime(), Binder);
801 if ((atom1->getType() != NULL)
802 && (atom1->getType()->getAtomicNumber() != 1)
803 && (atom2->getType() != NULL)
804 && (atom2->getType()->getAtomicNumber() != 1))
805 NoNonBonds++;
806
807 return Binder;
808};
809
810/** Set molecule::name from the basename without suffix in the given \a *filename.
811 * \param *filename filename
812 */
813void molecule::SetNameFromFilename(const char *filename)
814{
815 OBSERVE;
816 NOTIFY(MoleculeNameChanged);
817 int length = 0;
818 const char *molname = strrchr(filename, '/');
819 if (molname != NULL)
820 molname += sizeof(char); // search for filename without dirs
821 else
822 molname = filename; // contains no slashes
823 const char *endname = strchr(molname, '.');
824 if ((endname == NULL) || (endname < molname))
825 length = strlen(molname);
826 else
827 length = strlen(molname) - strlen(endname);
828 cout << "Set name of molecule " << getId() << " to " << molname << endl;
829 strncpy(name, molname, length);
830 name[length]='\0';
831};
832
833/** Removes atom from molecule list, but does not delete it.
834 * \param *pointer atom to be removed
835 * \return true - succeeded, false - atom not found in list
836 */
837bool molecule::UnlinkAtom(atom *pointer)
838{
839 if (pointer == NULL)
840 return false;
841 pointer->removeFromMolecule();
842 return true;
843};
844
845/** Removes every atom from molecule list.
846 * \return true - succeeded, false - atom not found in list
847 */
848bool molecule::CleanupMolecule()
849{
850 for (molecule::iterator iter = begin(); !empty(); iter = begin())
851 (*iter)->removeFromMolecule();
852 return empty();
853};
854
855/** Finds an atom specified by its continuous number.
856 * \param Nr number of atom withim molecule
857 * \return pointer to atom or NULL
858 */
859atom * molecule::FindAtom(int Nr) const
860{
861 LocalToGlobalId_t::const_iterator iter = LocalToGlobalId.find(Nr);
862 if (iter != LocalToGlobalId.end()) {
863 //LOG(0, "Found Atom Nr. " << walker->getNr());
864 return iter->second;
865 } else {
866 ELOG(1, "Atom with Nr " << Nr << " not found in molecule " << getName() << "'s list.");
867 return NULL;
868 }
869}
870
871/** Checks whether the given atom is a member of this molecule.
872 *
873 * We make use here of molecule::atomIds to get a result on
874 *
875 * @param _atom atom to check
876 * @return true - is member, false - is not
877 */
878bool molecule::isInMolecule(const atom * const _atom) const
879{
880 ASSERT(_atom->getMolecule() == this,
881 "molecule::isInMolecule() - atom is not designated to be in molecule '"
882 +toString(this->getName())+"'.");
883 molecule::const_iterator iter = atomIds.find(_atom->getId());
884 return (iter != atomIds.end());
885}
886
887/** Asks for atom number, and checks whether in list.
888 * \param *text question before entering
889 */
890atom * molecule::AskAtom(std::string text)
891{
892 int No;
893 atom *ion = NULL;
894 do {
895 //std::cout << "============Atom list==========================" << std::endl;
896 //mol->Output((ofstream *)&cout);
897 //std::cout << "===============================================" << std::endl;
898 std::cout << text;
899 cin >> No;
900 ion = this->FindAtom(No);
901 } while (ion == NULL);
902 return ion;
903};
904
905/** Checks if given coordinates are within cell volume.
906 * \param *x array of coordinates
907 * \return true - is within, false - out of cell
908 */
909bool molecule::CheckBounds(const Vector *x) const
910{
911 const RealSpaceMatrix &domain = World::getInstance().getDomain().getM();
912 bool result = true;
913 for (int i=0;i<NDIM;i++) {
914 result = result && ((x->at(i) >= 0) && (x->at(i) < domain.at(i,i)));
915 }
916 //return result;
917 return true; /// probably not gonna use the check no more
918};
919
920/** Prints molecule to *out.
921 * \param *out output stream
922 */
923bool molecule::Output(ostream * const output) const
924{
925 if (output == NULL) {
926 return false;
927 } else {
928 int AtomNo[MAX_ELEMENTS];
929 memset(AtomNo,0,(MAX_ELEMENTS-1)*sizeof(*AtomNo));
930 enumeration<const element*> elementLookup = formula.enumerateElements();
931 *output << "#Ion_TypeNr._Nr.R[0] R[1] R[2] MoveType (0 MoveIon, 1 FixedIon)" << endl;
932 for_each(begin(),end(),boost::bind(&atom::OutputArrayIndexed,_1,output,elementLookup,AtomNo,(const char*)0));
933 return true;
934 }
935};
936
937/** Outputs contents of each atom::ListOfBonds.
938 * \param *out output stream
939 */
940void molecule::OutputListOfBonds() const
941{
942 std::stringstream output;
943 LOG(2, "From Contents of ListOfBonds, all atoms:");
944 for (molecule::const_iterator iter = begin();
945 iter != end();
946 ++iter) {
947 (*iter)->OutputBondOfAtom(output);
948 output << std::endl << "\t\t";
949 }
950 LOG(2, output.str());
951}
952
953/** Returns an index map for two father-son-molecules.
954 * The map tells which atom in this molecule corresponds to which one in the other molecul with their fathers.
955 * \param *out output stream for debugging
956 * \param *OtherMolecule corresponding molecule with fathers
957 * \return allocated map of size molecule::AtomCount with map
958 * \todo make this with a good sort O(n), not O(n^2)
959 */
960int * molecule::GetFatherSonAtomicMap(const molecule * const OtherMolecule)
961{
962 LOG(3, "Begin of GetFatherAtomicMap.");
963 int *AtomicMap = new int[getAtomCount()];
964 for (int i=getAtomCount();i--;)
965 AtomicMap[i] = -1;
966 if (OtherMolecule == this) { // same molecule
967 for (int i=getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence
968 AtomicMap[i] = i;
969 LOG(4, "Map is trivial.");
970 } else {
971 std::stringstream output;
972 output << "Map is ";
973 for (molecule::const_iterator iter = const_cast<const molecule &>(*this).begin();
974 iter != const_cast<const molecule &>(*this).end();
975 ++iter) {
976 if ((*iter)->getFather() == NULL) {
977 AtomicMap[(*iter)->getNr()] = -2;
978 } else {
979 for (molecule::const_iterator runner = OtherMolecule->begin(); runner != OtherMolecule->end(); ++runner) {
980 //for (int i=0;i<AtomCount;i++) { // search atom
981 //for (int j=0;j<OtherMolecule->getAtomCount();j++) {
982 //LOG(4, "Comparing father " << (*iter)->getFather() << " with the other one " << (*runner)->getFather() << ".");
983 if ((*iter)->getFather() == (*runner))
984 AtomicMap[(*iter)->getNr()] = (*runner)->getNr();
985 }
986 }
987 output << AtomicMap[(*iter)->getNr()] << "\t";
988 }
989 LOG(4, output.str());
990 }
991 LOG(3, "End of GetFatherAtomicMap.");
992 return AtomicMap;
993};
994
995
996void molecule::flipActiveFlag(){
997 ActiveFlag = !ActiveFlag;
998}
999
1000Shape molecule::getBoundingShape(const double scale) const
1001{
1002 // create Sphere around every atom
1003 if (empty())
1004 return Nowhere();
1005 const_iterator iter = begin();
1006 const Vector center = (*iter)->getPosition();
1007 const double vdWRadius = (*iter)->getElement().getVanDerWaalsRadius();
1008 Shape BoundingShape = Sphere(center, vdWRadius*scale);
1009 for(++iter; iter != end(); ++iter) {
1010 const Vector center = (*iter)->getPosition();
1011 const double vdWRadius = (*iter)->getElement().getVanDerWaalsRadius();
1012 if (vdWRadius*scale != 0.)
1013 BoundingShape = Sphere(center, vdWRadius*scale) || BoundingShape;
1014 }
1015 return BoundingShape;
1016}
1017
1018Shape molecule::getBoundingSphere(const double boundary) const
1019{
1020 // get center and radius
1021 Vector center;
1022 double radius = 0.;
1023 {
1024 center.Zero();
1025 for(const_iterator iter = begin(); iter != end(); ++iter)
1026 center += (*iter)->getPosition();
1027 if (begin() != end())
1028 center *= 1./(double)size();
1029 for(const_iterator iter = begin(); iter != end(); ++iter) {
1030 const Vector &position = (*iter)->getPosition();
1031 const double temp_distance = position.DistanceSquared(center);
1032 if (temp_distance > radius)
1033 radius = temp_distance;
1034 }
1035 }
1036 // convert radius to true value and add some small boundary
1037 radius = sqrt(radius) + boundary + 1e+6*std::numeric_limits<double>::epsilon();
1038 LOG(1, "INFO: The " << size() << " atoms of the molecule are contained in a sphere at "
1039 << center << " with radius " << radius << ".");
1040
1041 // TODO: When we do not use a Sphere here anymore, then FillRegularGridAction will
1042 // will not work as it expects a sphere due to possible random rotations.
1043 Shape BoundingShape(Sphere(center, radius));
1044 LOG(1, "INFO: Created sphere at " << BoundingShape.getCenter() << " and radius "
1045 << BoundingShape.getRadius() << ".");
1046 return BoundingShape;
1047}
1048
1049molecule::BoundingBoxInfo molecule::updateBoundingBox() const
1050{
1051 BoundingBoxInfo info;
1052 Vector min = zeroVec;
1053 Vector max = zeroVec;
1054 for (int i=0;i<NDIM;++i) {
1055 if (!BoundingBoxSweepingAxis[i].right.empty()) {
1056 min[i] = BoundingBoxSweepingAxis[i].right.begin()->first;
1057 max[i] = BoundingBoxSweepingAxis[i].right.rbegin()->first;
1058 }
1059 }
1060 info.radius = (.5*(max-min)).Norm();
1061 info.position = .5*(max+min);
1062 return info;
1063}
1064
1065Vector molecule::updateMoleculeCenter() const
1066{
1067 return (1./(double)getAtomCount())*molcenter;
1068}
1069
1070molecule::BoundingBoxInfo molecule::getBoundingBox() const
1071{
1072 return **BoundingBox;
1073}
1074
1075Vector molecule::getMoleculeCenter() const
1076{
1077 return **MoleculeCenter;
1078}
1079
1080void molecule::update(Observable *publisher)
1081{
1082 ASSERT(0, "molecule::update() - did not sign on for any general updates.");
1083}
1084
1085void molecule::recieveNotification(Observable *publisher, Notification_ptr notification)
1086{
1087 const atom * const _atom = dynamic_cast<atom *>(publisher);
1088 if ((_atom != NULL) && containsAtom(_atom)) {
1089#ifdef LOG_OBSERVER
1090 observerLog().addMessage() << "++ Update of Observer "<< observerLog().getName(static_cast<Observer *>(this))
1091 << " received notification from atom " << _atom->getId() << " for channel "
1092 << notification->getChannelNo() << ".";
1093#endif
1094 switch (notification->getChannelNo()) {
1095 case AtomObservable::PositionChanged:
1096 {
1097 // emit others about one of our atoms moved
1098 _lastchangedatomid = _atom->getId();
1099 // update entry in map and also molecule center
1100 BoundingBoxInfo oldinfo = updateBoundingBox();
1101 for (int i=0;i<NDIM;++i) {
1102 AtomDistanceMap_t::left_iterator iter = BoundingBoxSweepingAxis[i].left.find(_atom->getId());
1103 ASSERT(iter != BoundingBoxSweepingAxis[i].left.end(),
1104 "molecule::recieveNotification() - could not find atom "+toString(_atom->getId())
1105 +" in BoundingBoxSweepingAxis.");
1106 molcenter[i] -= iter->second;
1107 BoundingBoxSweepingAxis[i].left.erase(iter);
1108 const Vector &position = _atom->getPosition();
1109 BoundingBoxSweepingAxis[i].left.insert(
1110 std::make_pair(_atom->getId(), position[i]) );
1111 molcenter[i] += position[i];
1112 }
1113 BoundingBoxInfo newinfo = updateBoundingBox();
1114 OBSERVE;
1115 NOTIFY(AtomMoved);
1116 NOTIFY(MoleculeCenterChanged);
1117 if (oldinfo != newinfo)
1118 NOTIFY(BoundingBoxChanged);
1119 break;
1120 }
1121 case AtomObservable::ElementChanged:
1122 {
1123 // emit others about one of our atoms moved
1124 _lastchangedatomid = _atom->getId();
1125 OBSERVE;
1126 NOTIFY(FormulaChanged);
1127 const ElementPerAtom_t::iterator iter = ElementPerAtom.find(_lastchangedatomid);
1128 ASSERT( iter != ElementPerAtom.end(),
1129 "molecule::recieveNotification() - atom "
1130 +toString(_atom->getId()+" is not contained in ElementsPerAtom."));
1131 formula -= iter->second;
1132 if (iter->second == (atomicNumber_t)1) // was a hydrogen ?
1133 --NoNonHydrogen;
1134 iter->second = _atom->getElementNo();
1135 formula += iter->second;
1136 if (iter->second == (atomicNumber_t)1) // is a hydrogen ?
1137 ++NoNonHydrogen;
1138 break;
1139 }
1140 case AtomObservable::BondsAdded:
1141 case AtomObservable::BondsRemoved:
1142 {
1143 // emit others about one of our atoms moved
1144 _lastchangedatomid = _atom->getId();
1145 const BondCountsPerAtom_t::iterator iter = BondCountsPerAtom.find(_lastchangedatomid);
1146 ASSERT( iter != BondCountsPerAtom.end(),
1147 "molecule::recieveNotification() - atom "
1148 +toString(_atom->getId()+" is not contained in BondCountsPerAtom."));
1149 BondCount -= iter->second;
1150 iter->second = _atom->getListOfBonds().size();
1151 BondCount += iter->second;
1152 break;
1153 }
1154 default:
1155 ASSERT( 0, "molecule::recieveNotification() - we did not sign up for channel "
1156 +toString(notification->getChannelNo()));
1157 break;
1158 }
1159 }
1160}
1161
1162void molecule::subjectKilled(Observable *publisher)
1163{
1164 // do nothing, atom does it all
1165}
1166
1167void molecule::select()
1168{
1169 OBSERVE;
1170 selected = true;
1171 NOTIFY(SelectionChanged);
1172}
1173
1174void molecule::unselect()
1175{
1176 OBSERVE;
1177 selected = false;
1178 NOTIFY(SelectionChanged);
1179}
1180
1181void molecule::associateAtomWithMolecule(atom *_atom)
1182{
1183 _atom->signOn(this, AtomObservable::PositionChanged);
1184 _atom->signOn(this, AtomObservable::ElementChanged);
1185 _atom->signOn(this, AtomObservable::BondsAdded);
1186 _atom->signOn(this, AtomObservable::BondsRemoved);
1187 insert(_atom);
1188 {
1189 const size_t atom_bondcount = _atom->getListOfBonds().size();
1190#ifndef NDEBUG
1191 const std::pair<BondCountsPerAtom_t::iterator, bool> inserter =
1192#endif
1193 BondCountsPerAtom.insert( std::make_pair(_atom->getId(), atom_bondcount) );
1194 ASSERT( inserter.second,
1195 "molecule::associateAtomWithMolecule() - atom "
1196 +toString(_atom->getId()+" already in BondCountsPerAtom."));
1197 BondCount += atom_bondcount;
1198 }
1199 {
1200 OBSERVE;
1201 const int atom_elementno = _atom->getElementNo();
1202#ifndef NDEBUG
1203 const std::pair<ElementPerAtom_t::iterator, bool> inserter =
1204#endif
1205 ElementPerAtom.insert( std::make_pair(_atom->getId(), atom_elementno) );
1206 ASSERT( inserter.second,
1207 "molecule::associateAtomWithMolecule() - atom "
1208 +toString(_atom->getId()+" already in ElementPerAtom."));
1209 formula += atom_elementno;
1210 NOTIFY(FormulaChanged);
1211 }
1212}
1213
1214void molecule::disassociateAtomWithMolecule(atom *_atom)
1215{
1216 _atom->signOff(this, AtomObservable::PositionChanged);
1217 _atom->signOff(this, AtomObservable::ElementChanged);
1218 _atom->signOff(this, AtomObservable::BondsAdded);
1219 _atom->signOff(this, AtomObservable::BondsRemoved);
1220 erase(_atom);
1221 {
1222 const BondCountsPerAtom_t::iterator iter = BondCountsPerAtom.find(_atom->getId());
1223 ASSERT( iter != BondCountsPerAtom.end(),
1224 "molecule::disassociateAtomWithMolecule() - atom "
1225 +toString(_atom->getId()+" is not contained in BondCountsPerAtom."));
1226 BondCount -= iter->second;
1227 BondCountsPerAtom.erase(iter);
1228 }
1229 {
1230 OBSERVE;
1231 const ElementPerAtom_t::iterator iter = ElementPerAtom.find(_atom->getId());
1232 ASSERT( iter != ElementPerAtom.end(),
1233 "molecule::disassociateAtomWithMolecule() - atom "
1234 +toString(_atom->getId()+" is not contained in ElementPerAtom."));
1235 formula -= iter->second;
1236 ElementPerAtom.erase(iter);
1237 NOTIFY(FormulaChanged);
1238 }
1239}
1240
1241// construct idpool
1242CONSTRUCT_IDPOOL(atomId_t, continuousId)
1243
Note: See TracBrowser for help on using the repository browser.