[bcf653] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
[0aa122] | 4 | * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
|
---|
[5aaa43] | 5 | * Copyright (C) 2013 Frederik Heber. All rights reserved.
|
---|
[94d5ac6] | 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/>.
|
---|
[bcf653] | 22 | */
|
---|
| 23 |
|
---|
[14de469] | 24 | /** \file molecules.cpp
|
---|
[69eb71] | 25 | *
|
---|
[14de469] | 26 | * Functions for the class molecule.
|
---|
[69eb71] | 27 | *
|
---|
[14de469] | 28 | */
|
---|
| 29 |
|
---|
[bf3817] | 30 | // include config.h
|
---|
[aafd77] | 31 | #ifdef HAVE_CONFIG_H
|
---|
| 32 | #include <config.h>
|
---|
| 33 | #endif
|
---|
| 34 |
|
---|
[ad011c] | 35 | #include "CodePatterns/MemDebug.hpp"
|
---|
[112b09] | 36 |
|
---|
[0a5beb] | 37 | #include <algorithm>
|
---|
[e39e7a] | 38 | #include <boost/assign.hpp>
|
---|
[ac9b56] | 39 | #include <boost/bind.hpp>
|
---|
[9df5c6] | 40 | #include <boost/foreach.hpp>
|
---|
[0a5beb] | 41 | #include <cstring>
|
---|
[49e1ae] | 42 |
|
---|
[aafd77] | 43 | #include <gsl/gsl_inline.h>
|
---|
| 44 | #include <gsl/gsl_heapsort.h>
|
---|
| 45 |
|
---|
[560bbe] | 46 | #include "molecule.hpp"
|
---|
| 47 |
|
---|
[6f0841] | 48 | #include "Atom/atom.hpp"
|
---|
[129204] | 49 | #include "Bond/bond.hpp"
|
---|
[9d83b6] | 50 | #include "Box.hpp"
|
---|
| 51 | #include "CodePatterns/enumeration.hpp"
|
---|
| 52 | #include "CodePatterns/Log.hpp"
|
---|
[e39e7a] | 53 | #include "CodePatterns/Observer/Observable.hpp"
|
---|
[c32d21] | 54 | #include "CodePatterns/Observer/Notification.hpp"
|
---|
[a80fbdf] | 55 | #include "config.hpp"
|
---|
[560bbe] | 56 | #include "Descriptors/AtomIdDescriptor.hpp"
|
---|
[3bdb6d] | 57 | #include "Element/element.hpp"
|
---|
[129204] | 58 | #include "Graph/BondGraph.hpp"
|
---|
[783e88] | 59 | #include "LinearAlgebra/Exceptions.hpp"
|
---|
[13d150] | 60 | #include "LinearAlgebra/leastsquaremin.hpp"
|
---|
[9d83b6] | 61 | #include "LinearAlgebra/Plane.hpp"
|
---|
| 62 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
| 63 | #include "LinearAlgebra/Vector.hpp"
|
---|
[53c7fc] | 64 | #include "LinkedCell/linkedcell.hpp"
|
---|
[5d8f4f] | 65 | #include "MoleculeObserver.hpp"
|
---|
[560bbe] | 66 | #include "IdPool_impl.hpp"
|
---|
[c67ff9] | 67 | #include "Shapes/BaseShapes.hpp"
|
---|
[d127c8] | 68 | #include "Tesselation/tesselation.hpp"
|
---|
[b34306] | 69 | #include "World.hpp"
|
---|
[9d83b6] | 70 | #include "WorldTime.hpp"
|
---|
[14de469] | 71 |
|
---|
[e39e7a] | 72 | using namespace boost::assign;
|
---|
| 73 |
|
---|
| 74 | // static entities
|
---|
[cbd409] | 75 | static Observable::channels_t getAtomPositionsChannels()
|
---|
[e39e7a] | 76 | {
|
---|
| 77 | Observable::channels_t channels;
|
---|
| 78 | channels += molecule::AtomInserted, molecule::AtomRemoved, molecule::AtomMoved;
|
---|
| 79 | return channels;
|
---|
| 80 | }
|
---|
[14de469] | 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 | */
|
---|
[4d2b33] | 87 | molecule::molecule() :
|
---|
[cd5047] | 88 | Observable("molecule"),
|
---|
[458c31] | 89 | MDSteps(0),
|
---|
| 90 | NoNonBonds(0),
|
---|
| 91 | NoCyclicBonds(0),
|
---|
| 92 | ActiveFlag(false),
|
---|
| 93 | IndexNr(-1),
|
---|
[29f7c1] | 94 | NoNonHydrogen(0),
|
---|
| 95 | BondCount(0),
|
---|
[52ed5b] | 96 | atomIdPool(1, 20, 100),
|
---|
[e39e7a] | 97 | BoundingBoxSweepingAxis(std::vector<AtomDistanceMap_t>(NDIM)),
|
---|
[fb95a5] | 98 | _lastchangedatomid(-1),
|
---|
[cbd409] | 99 | last_atom(0),
|
---|
[b71881] | 100 | molcenter(zeroVec),
|
---|
| 101 | selected(false)
|
---|
[69eb71] | 102 | {
|
---|
[6a3c83] | 103 | // add specific channels
|
---|
| 104 | Channels *OurChannel = new Channels;
|
---|
[574d377] | 105 | Observable::insertNotificationChannel( std::make_pair( static_cast<Observable *>(this), OurChannel) );
|
---|
[6a3c83] | 106 | for (size_t type = 0; type < (size_t)NotificationType_MAX; ++type)
|
---|
| 107 | OurChannel->addChannel(type);
|
---|
[fa649a] | 108 |
|
---|
[e39e7a] | 109 | // cannot initialize in initializer body as then channels have not been setup yet
|
---|
| 110 | BoundingBox.reset(
|
---|
| 111 | new Cacheable<BoundingBoxInfo>(
|
---|
[cbd409] | 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()));
|
---|
[e39e7a] | 116 |
|
---|
[387b36] | 117 | strcpy(name,World::getInstance().getDefaultName().c_str());
|
---|
[5d8f4f] | 118 |
|
---|
| 119 | // inform MoleculeObserver about new molecule
|
---|
[d4ba3f] | 120 | MoleculeObserver::getInstance().Inserted(this);
|
---|
[e39e7a] | 121 | }
|
---|
[14de469] | 122 |
|
---|
[cbc5fb] | 123 | molecule *NewMolecule(){
|
---|
[4d2b33] | 124 | return new molecule();
|
---|
[cbc5fb] | 125 | }
|
---|
| 126 |
|
---|
[14de469] | 127 | /** Destructor of class molecule.
|
---|
| 128 | * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
|
---|
| 129 | */
|
---|
[69eb71] | 130 | molecule::~molecule()
|
---|
[14de469] | 131 | {
|
---|
[042f82] | 132 | CleanupMolecule();
|
---|
[5d8f4f] | 133 |
|
---|
| 134 | // inform MoleculeObserver about removed molecule
|
---|
[d4ba3f] | 135 | MoleculeObserver::getInstance().Removed(this);
|
---|
[14de469] | 136 | };
|
---|
| 137 |
|
---|
[357fba] | 138 |
|
---|
[cbc5fb] | 139 | void DeleteMolecule(molecule *mol){
|
---|
| 140 | delete mol;
|
---|
| 141 | }
|
---|
| 142 |
|
---|
[520c8b] | 143 | // getter and setter
|
---|
[73a857] | 144 | const std::string molecule::getName() const{
|
---|
[520c8b] | 145 | return std::string(name);
|
---|
| 146 | }
|
---|
| 147 |
|
---|
[ea7176] | 148 | int molecule::getAtomCount() const{
|
---|
[e791dc] | 149 | return atomIds.size();
|
---|
| 150 | }
|
---|
| 151 |
|
---|
[520c8b] | 152 | void molecule::setName(const std::string _name){
|
---|
[2ba827] | 153 | OBSERVE;
|
---|
[6a3c83] | 154 | NOTIFY(MoleculeNameChanged);
|
---|
[35b698] | 155 | cout << "Set name of molecule " << getId() << " to " << _name << endl;
|
---|
[520c8b] | 156 | strncpy(name,_name.c_str(),MAXSTRINGSIZE);
|
---|
| 157 | }
|
---|
| 158 |
|
---|
[c6ab91] | 159 | void 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 |
|
---|
[560bbe] | 169 | bool molecule::changeAtomNr(int oldNr, int newNr, atom* target){
|
---|
| 170 | OBSERVE;
|
---|
| 171 | if(atomIdPool.reserveId(newNr)){
|
---|
[6a3c83] | 172 | NOTIFY(AtomNrChanged);
|
---|
[560bbe] | 173 | if (oldNr != -1) // -1 is reserved and indicates no number
|
---|
| 174 | atomIdPool.releaseId(oldNr);
|
---|
[c6ab91] | 175 | LocalToGlobalId.erase(oldNr);
|
---|
[560bbe] | 176 | ASSERT (target,
|
---|
| 177 | "molecule::changeAtomNr() - given target is NULL, cannot set Nr or name.");
|
---|
| 178 | target->setNr(newNr);
|
---|
[fb95a5] | 179 | _lastchangedatomid = target->getId();
|
---|
[c6ab91] | 180 | InsertLocalToGlobalId(target);
|
---|
[560bbe] | 181 | setAtomName(target);
|
---|
| 182 | return true;
|
---|
| 183 | } else{
|
---|
| 184 | return false;
|
---|
| 185 | }
|
---|
| 186 | }
|
---|
| 187 |
|
---|
[ceaab1] | 188 | bool 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 );
|
---|
[e39e7a] | 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 | }
|
---|
[ceaab1] | 205 | return true;
|
---|
| 206 | }
|
---|
| 207 |
|
---|
[a7a087] | 208 | bool 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)){
|
---|
[f54524] | 212 | OBSERVE;
|
---|
| 213 | NOTIFY(IndexChanged);
|
---|
[a7a087] | 214 | id = newId;
|
---|
| 215 | return true;
|
---|
| 216 | }
|
---|
| 217 | else{
|
---|
| 218 | return false;
|
---|
| 219 | }
|
---|
| 220 | }
|
---|
| 221 |
|
---|
| 222 |
|
---|
[73a857] | 223 | moleculeId_t molecule::getId() const {
|
---|
[cbc5fb] | 224 | return id;
|
---|
| 225 | }
|
---|
| 226 |
|
---|
| 227 | void molecule::setId(moleculeId_t _id){
|
---|
| 228 | id =_id;
|
---|
| 229 | }
|
---|
| 230 |
|
---|
[73a857] | 231 | const Formula &molecule::getFormula() const {
|
---|
[f17e1c] | 232 | return formula;
|
---|
[ac9b56] | 233 | }
|
---|
| 234 |
|
---|
[73a857] | 235 | unsigned int molecule::getElementCount() const{
|
---|
[389cc8] | 236 | return formula.getElementCount();
|
---|
| 237 | }
|
---|
| 238 |
|
---|
| 239 | bool molecule::hasElement(const element *element) const{
|
---|
| 240 | return formula.hasElement(element);
|
---|
| 241 | }
|
---|
| 242 |
|
---|
| 243 | bool molecule::hasElement(atomicNumber_t Z) const{
|
---|
| 244 | return formula.hasElement(Z);
|
---|
| 245 | }
|
---|
| 246 |
|
---|
| 247 | bool molecule::hasElement(const string &shorthand) const{
|
---|
| 248 | return formula.hasElement(shorthand);
|
---|
| 249 | }
|
---|
| 250 |
|
---|
[bd58fb] | 251 | /************************** Access to the List of Atoms ****************/
|
---|
| 252 |
|
---|
[9879f6] | 253 | molecule::const_iterator molecule::erase( const_iterator loc )
|
---|
| 254 | {
|
---|
[bf8e20] | 255 | OBSERVE;
|
---|
[59fff1] | 256 | const_iterator iter = loc;
|
---|
[30c753] | 257 | ++iter;
|
---|
[59fff1] | 258 | atom * const _atom = const_cast<atom *>(*loc);
|
---|
[8c001a] | 259 | {
|
---|
[fb95a5] | 260 | _lastchangedatomid = _atom->getId();
|
---|
[8c001a] | 261 | NOTIFY(AtomRemoved);
|
---|
| 262 | }
|
---|
[59fff1] | 263 | atomIds.erase( _atom->getId() );
|
---|
[e39e7a] | 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 | }
|
---|
[cbd409] | 272 | {
|
---|
| 273 | molcenter -= _atom->getPosition();
|
---|
| 274 | }
|
---|
[6a3c83] | 275 | {
|
---|
| 276 | NOTIFY(AtomNrChanged);
|
---|
| 277 | atomIdPool.releaseId(_atom->getNr());
|
---|
[c6ab91] | 278 | LocalToGlobalId.erase(_atom->getNr());
|
---|
[6a3c83] | 279 | _atom->setNr(-1);
|
---|
| 280 | }
|
---|
[59fff1] | 281 | _atom->removeFromMolecule();
|
---|
[9879f6] | 282 | return iter;
|
---|
| 283 | }
|
---|
| 284 |
|
---|
[6cfa36] | 285 | molecule::const_iterator molecule::erase( atom * key )
|
---|
[9879f6] | 286 | {
|
---|
[f01769] | 287 | const_iterator iter = const_cast<const molecule &>(*this).find(key);
|
---|
[a063787] | 288 | if (iter != const_cast<const molecule &>(*this).end())
|
---|
| 289 | return erase(iter);
|
---|
| 290 | else
|
---|
| 291 | return iter;
|
---|
[9879f6] | 292 | }
|
---|
| 293 |
|
---|
| 294 | pair<molecule::iterator,bool> molecule::insert ( atom * const key )
|
---|
| 295 | {
|
---|
[bf8e20] | 296 | OBSERVE;
|
---|
[6a3c83] | 297 | NOTIFY(AtomInserted);
|
---|
[fb95a5] | 298 | _lastchangedatomid = key->getId();
|
---|
[8e1f901] | 299 | std::pair<iterator,bool> res = atomIds.insert(key->getId());
|
---|
[274d45] | 300 | if (res.second) { // push atom if went well
|
---|
[e39e7a] | 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 | }
|
---|
[cbd409] | 309 | {
|
---|
| 310 | molcenter += key->getPosition();
|
---|
| 311 | }
|
---|
[6a3c83] | 312 | NOTIFY(AtomNrChanged);
|
---|
[560bbe] | 313 | key->setNr(atomIdPool.getNextId());
|
---|
[c6ab91] | 314 | InsertLocalToGlobalId(key);
|
---|
[560bbe] | 315 | setAtomName(key);
|
---|
[8e1f901] | 316 | return res;
|
---|
[274d45] | 317 | } else {
|
---|
[30c753] | 318 | return pair<iterator,bool>(end(),res.second);
|
---|
[274d45] | 319 | }
|
---|
[9879f6] | 320 | }
|
---|
[520c8b] | 321 |
|
---|
[560bbe] | 322 | void molecule::setAtomName(atom *_atom) const
|
---|
| 323 | {
|
---|
| 324 | std::stringstream sstr;
|
---|
[52ed5b] | 325 | sstr << _atom->getType()->getSymbol() << _atom->getNr();
|
---|
[560bbe] | 326 | _atom->setName(sstr.str());
|
---|
| 327 | }
|
---|
| 328 |
|
---|
[f01769] | 329 | World::AtomComposite molecule::getAtomSet()
|
---|
[3738f0] | 330 | {
|
---|
[9317be] | 331 | World::AtomComposite vector_of_atoms;
|
---|
[59fff1] | 332 | for (molecule::iterator iter = begin(); iter != end(); ++iter)
|
---|
[30c753] | 333 | vector_of_atoms.push_back(*iter);
|
---|
[3738f0] | 334 | return vector_of_atoms;
|
---|
| 335 | }
|
---|
| 336 |
|
---|
[f01769] | 337 | World::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 |
|
---|
[14de469] | 345 | /** Adds given atom \a *pointer from molecule list.
|
---|
[69eb71] | 346 | * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount
|
---|
[14de469] | 347 | * \param *pointer allocated and set atom
|
---|
| 348 | * \return true - succeeded, false - atom not found in list
|
---|
| 349 | */
|
---|
| 350 | bool molecule::AddAtom(atom *pointer)
|
---|
[69eb71] | 351 | {
|
---|
[042f82] | 352 | if (pointer != NULL) {
|
---|
[356ae4] | 353 | // molecule::insert() is called by setMolecule()
|
---|
[6cfa36] | 354 | pointer->setMolecule(this);
|
---|
[f721c6] | 355 | }
|
---|
[9879f6] | 356 | return true;
|
---|
[14de469] | 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
|
---|
[89c8b2] | 362 | * \return pointer to the newly added atom
|
---|
[14de469] | 363 | */
|
---|
| 364 | atom * molecule::AddCopyAtom(atom *pointer)
|
---|
[69eb71] | 365 | {
|
---|
[f721c6] | 366 | atom *retval = NULL;
|
---|
[042f82] | 367 | if (pointer != NULL) {
|
---|
[46d958] | 368 | atom *walker = pointer->clone();
|
---|
[c6ab91] | 369 | AddAtom(walker);
|
---|
[f721c6] | 370 | retval=walker;
|
---|
| 371 | }
|
---|
| 372 | return retval;
|
---|
[14de469] | 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
|
---|
[042f82] | 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.
|
---|
[14de469] | 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
|
---|
[042f82] | 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$.
|
---|
[69eb71] | 399 | *
|
---|
[14de469] | 400 | * \param *out output stream for debugging
|
---|
[69eb71] | 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)
|
---|
[14de469] | 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 | */
|
---|
[06804b] | 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 | //};
|
---|
[14de469] | 656 |
|
---|
| 657 | /** Creates a copy of this molecule.
|
---|
[c67ff9] | 658 | * \param offset translation Vector for the new molecule relative to old one
|
---|
[14de469] | 659 | * \return copy of molecule
|
---|
| 660 | */
|
---|
[f01769] | 661 | molecule *molecule::CopyMolecule(const Vector &offset)
|
---|
[14de469] | 662 | {
|
---|
[5f612ee] | 663 | molecule *copy = World::getInstance().createMolecule();
|
---|
[042f82] | 664 |
|
---|
| 665 | // copy all atoms
|
---|
[30c753] | 666 | std::map< const atom *, atom *> FatherFinder;
|
---|
[59fff1] | 667 | for (iterator iter = begin(); iter != end(); ++iter) {
|
---|
| 668 | atom * const copy_atom = copy->AddCopyAtom(*iter);
|
---|
[c67ff9] | 669 | copy_atom->setPosition(copy_atom->getPosition() + offset);
|
---|
[30c753] | 670 | FatherFinder.insert( std::make_pair( *iter, copy_atom ) );
|
---|
| 671 | }
|
---|
[042f82] | 672 |
|
---|
| 673 | // copy all bonds
|
---|
[f01769] | 674 | for(const_iterator AtomRunner = const_cast<const molecule &>(*this).begin();
|
---|
| 675 | AtomRunner != const_cast<const molecule &>(*this).end();
|
---|
| 676 | ++AtomRunner) {
|
---|
[9d83b6] | 677 | const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
|
---|
| 678 | for(BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
| 679 | BondRunner != ListOfBonds.end();
|
---|
| 680 | ++BondRunner)
|
---|
[e08c46] | 681 | if ((*BondRunner)->leftatom == *AtomRunner) {
|
---|
[88c8ec] | 682 | bond::ptr Binder = (*BondRunner);
|
---|
[e08c46] | 683 | // get the pendant atoms of current bond in the copy molecule
|
---|
[30c753] | 684 | ASSERT(FatherFinder.count(Binder->leftatom),
|
---|
[59fff1] | 685 | "molecule::CopyMolecule() - No copy of original left atom "
|
---|
| 686 | +toString(Binder->leftatom)+" for bond copy found");
|
---|
[30c753] | 687 | ASSERT(FatherFinder.count(Binder->rightatom),
|
---|
[59fff1] | 688 | "molecule::CopyMolecule() - No copy of original right atom "
|
---|
| 689 | +toString(Binder->rightatom)+" for bond copy found");
|
---|
[30c753] | 690 | atom * const LeftAtom = FatherFinder[Binder->leftatom];
|
---|
| 691 | atom * const RightAtom = FatherFinder[Binder->rightatom];
|
---|
| 692 |
|
---|
[1f693d] | 693 | bond::ptr const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->getDegree());
|
---|
[e08c46] | 694 | NewBond->Cyclic = Binder->Cyclic;
|
---|
| 695 | if (Binder->Cyclic)
|
---|
| 696 | copy->NoCyclicBonds++;
|
---|
| 697 | NewBond->Type = Binder->Type;
|
---|
| 698 | }
|
---|
[9d83b6] | 699 | }
|
---|
[042f82] | 700 | // correct fathers
|
---|
[30c753] | 701 | //for_each(begin(),end(),mem_fun(&atom::CorrectFather));
|
---|
[cee0b57] | 702 |
|
---|
[042f82] | 703 | return copy;
|
---|
[14de469] | 704 | };
|
---|
| 705 |
|
---|
[89c8b2] | 706 |
|
---|
[9df680] | 707 | /** Destroys all atoms inside this molecule.
|
---|
| 708 | */
|
---|
[a7aebd] | 709 | void removeAtomsinMolecule(molecule *&_mol)
|
---|
[9df680] | 710 | {
|
---|
[0a5beb] | 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());
|
---|
[9df680] | 715 | // remove each atom from world
|
---|
[0a5beb] | 716 | for(std::vector<atom *>::iterator AtomRunner = atoms.begin();
|
---|
| 717 | AtomRunner != atoms.end(); ++AtomRunner)
|
---|
[9df680] | 718 | World::getInstance().destroyAtom(*AtomRunner);
|
---|
[a7aebd] | 719 | // make sure that pointer os not usable
|
---|
| 720 | _mol = NULL;
|
---|
[9df680] | 721 | };
|
---|
| 722 |
|
---|
| 723 |
|
---|
[89c8b2] | 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 | */
|
---|
[f01769] | 730 | molecule* molecule::CopyMoleculeFromSubRegion(const Shape ®ion) {
|
---|
[5f612ee] | 731 | molecule *copy = World::getInstance().createMolecule();
|
---|
[89c8b2] | 732 |
|
---|
[30c753] | 733 | // copy all atoms
|
---|
| 734 | std::map< const atom *, atom *> FatherFinder;
|
---|
[59fff1] | 735 | for (iterator iter = begin(); iter != end(); ++iter) {
|
---|
[4e904b] | 736 | if (region.isInside((*iter)->getPosition())) {
|
---|
[59fff1] | 737 | atom * const copy_atom = copy->AddCopyAtom(*iter);
|
---|
[30c753] | 738 | FatherFinder.insert( std::make_pair( *iter, copy_atom ) );
|
---|
[9df5c6] | 739 | }
|
---|
| 740 | }
|
---|
[89c8b2] | 741 |
|
---|
[30c753] | 742 | // copy all bonds
|
---|
[f01769] | 743 | for(molecule::const_iterator AtomRunner = const_cast<const molecule &>(*this).begin();
|
---|
| 744 | AtomRunner != const_cast<const molecule &>(*this).end();
|
---|
| 745 | ++AtomRunner) {
|
---|
[30c753] | 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) {
|
---|
[88c8ec] | 751 | bond::ptr Binder = (*BondRunner);
|
---|
[30c753] | 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 |
|
---|
[1f693d] | 758 | bond::ptr const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->getDegree());
|
---|
[30c753] | 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 |
|
---|
[e138de] | 769 | //TODO: copy->BuildInducedSubgraph(this);
|
---|
[89c8b2] | 770 |
|
---|
| 771 | return copy;
|
---|
| 772 | }
|
---|
| 773 |
|
---|
[14de469] | 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 | */
|
---|
[88c8ec] | 780 | bond::ptr molecule::AddBond(atom *atom1, atom *atom2, int degree)
|
---|
[14de469] | 781 | {
|
---|
[7d82a5] | 782 | bond::ptr Binder;
|
---|
[05a97c] | 783 |
|
---|
| 784 | // some checks to make sure we are able to create the bond
|
---|
[59fff1] | 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");
|
---|
[05a97c] | 797 |
|
---|
[7d82a5] | 798 | Binder.reset(new bond(atom1, atom2, degree));
|
---|
[073a9e4] | 799 | atom1->RegisterBond(WorldTime::getTime(), Binder);
|
---|
| 800 | atom2->RegisterBond(WorldTime::getTime(), Binder);
|
---|
[59fff1] | 801 | if ((atom1->getType() != NULL)
|
---|
| 802 | && (atom1->getType()->getAtomicNumber() != 1)
|
---|
| 803 | && (atom2->getType() != NULL)
|
---|
| 804 | && (atom2->getType()->getAtomicNumber() != 1))
|
---|
[05a97c] | 805 | NoNonBonds++;
|
---|
| 806 |
|
---|
[042f82] | 807 | return Binder;
|
---|
[14de469] | 808 | };
|
---|
| 809 |
|
---|
[1907a7] | 810 | /** Set molecule::name from the basename without suffix in the given \a *filename.
|
---|
| 811 | * \param *filename filename
|
---|
| 812 | */
|
---|
[d67150] | 813 | void molecule::SetNameFromFilename(const char *filename)
|
---|
[1907a7] | 814 | {
|
---|
[575343] | 815 | OBSERVE;
|
---|
[1c9eb8] | 816 | NOTIFY(MoleculeNameChanged);
|
---|
[1907a7] | 817 | int length = 0;
|
---|
[f7f7a4] | 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
|
---|
[49e1ae] | 823 | const char *endname = strchr(molname, '.');
|
---|
[1907a7] | 824 | if ((endname == NULL) || (endname < molname))
|
---|
| 825 | length = strlen(molname);
|
---|
| 826 | else
|
---|
| 827 | length = strlen(molname) - strlen(endname);
|
---|
[35b698] | 828 | cout << "Set name of molecule " << getId() << " to " << molname << endl;
|
---|
[1907a7] | 829 | strncpy(name, molname, length);
|
---|
[d67150] | 830 | name[length]='\0';
|
---|
[1907a7] | 831 | };
|
---|
| 832 |
|
---|
[cee0b57] | 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
|
---|
[f3278b] | 836 | */
|
---|
[cee0b57] | 837 | bool molecule::UnlinkAtom(atom *pointer)
|
---|
[f3278b] | 838 | {
|
---|
[cee0b57] | 839 | if (pointer == NULL)
|
---|
| 840 | return false;
|
---|
[2e4105] | 841 | pointer->removeFromMolecule();
|
---|
[cee0b57] | 842 | return true;
|
---|
[f3278b] | 843 | };
|
---|
| 844 |
|
---|
[cee0b57] | 845 | /** Removes every atom from molecule list.
|
---|
| 846 | * \return true - succeeded, false - atom not found in list
|
---|
[14de469] | 847 | */
|
---|
[cee0b57] | 848 | bool molecule::CleanupMolecule()
|
---|
[14de469] | 849 | {
|
---|
[9879f6] | 850 | for (molecule::iterator iter = begin(); !empty(); iter = begin())
|
---|
[2e4105] | 851 | (*iter)->removeFromMolecule();
|
---|
[274d45] | 852 | return empty();
|
---|
[69eb71] | 853 | };
|
---|
[14de469] | 854 |
|
---|
[cee0b57] | 855 | /** Finds an atom specified by its continuous number.
|
---|
| 856 | * \param Nr number of atom withim molecule
|
---|
| 857 | * \return pointer to atom or NULL
|
---|
[14de469] | 858 | */
|
---|
[9879f6] | 859 | atom * molecule::FindAtom(int Nr) const
|
---|
| 860 | {
|
---|
[c6ab91] | 861 | LocalToGlobalId_t::const_iterator iter = LocalToGlobalId.find(Nr);
|
---|
| 862 | if (iter != LocalToGlobalId.end()) {
|
---|
[47d041] | 863 | //LOG(0, "Found Atom Nr. " << walker->getNr());
|
---|
[c6ab91] | 864 | return iter->second;
|
---|
[cee0b57] | 865 | } else {
|
---|
[ca8bea] | 866 | ELOG(1, "Atom with Nr " << Nr << " not found in molecule " << getName() << "'s list.");
|
---|
[cee0b57] | 867 | return NULL;
|
---|
[042f82] | 868 | }
|
---|
[59fff1] | 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 | */
|
---|
[f01769] | 878 | bool molecule::isInMolecule(const atom * const _atom) const
|
---|
[59fff1] | 879 | {
|
---|
| 880 | ASSERT(_atom->getMolecule() == this,
|
---|
| 881 | "molecule::isInMolecule() - atom is not designated to be in molecule '"
|
---|
| 882 | +toString(this->getName())+"'.");
|
---|
[8e1f901] | 883 | molecule::const_iterator iter = atomIds.find(_atom->getId());
|
---|
[59fff1] | 884 | return (iter != atomIds.end());
|
---|
| 885 | }
|
---|
[14de469] | 886 |
|
---|
[cee0b57] | 887 | /** Asks for atom number, and checks whether in list.
|
---|
| 888 | * \param *text question before entering
|
---|
[a6b7fb] | 889 | */
|
---|
[955b91] | 890 | atom * molecule::AskAtom(std::string text)
|
---|
[a6b7fb] | 891 | {
|
---|
[cee0b57] | 892 | int No;
|
---|
| 893 | atom *ion = NULL;
|
---|
| 894 | do {
|
---|
[47d041] | 895 | //std::cout << "============Atom list==========================" << std::endl;
|
---|
[cee0b57] | 896 | //mol->Output((ofstream *)&cout);
|
---|
[47d041] | 897 | //std::cout << "===============================================" << std::endl;
|
---|
| 898 | std::cout << text;
|
---|
[cee0b57] | 899 | cin >> No;
|
---|
| 900 | ion = this->FindAtom(No);
|
---|
| 901 | } while (ion == NULL);
|
---|
| 902 | return ion;
|
---|
[a6b7fb] | 903 | };
|
---|
| 904 |
|
---|
[cee0b57] | 905 | /** Checks if given coordinates are within cell volume.
|
---|
| 906 | * \param *x array of coordinates
|
---|
| 907 | * \return true - is within, false - out of cell
|
---|
[14de469] | 908 | */
|
---|
[cee0b57] | 909 | bool molecule::CheckBounds(const Vector *x) const
|
---|
[14de469] | 910 | {
|
---|
[cca9ef] | 911 | const RealSpaceMatrix &domain = World::getInstance().getDomain().getM();
|
---|
[cee0b57] | 912 | bool result = true;
|
---|
| 913 | for (int i=0;i<NDIM;i++) {
|
---|
[84c494] | 914 | result = result && ((x->at(i) >= 0) && (x->at(i) < domain.at(i,i)));
|
---|
[042f82] | 915 | }
|
---|
[cee0b57] | 916 | //return result;
|
---|
| 917 | return true; /// probably not gonna use the check no more
|
---|
[69eb71] | 918 | };
|
---|
[14de469] | 919 |
|
---|
[cee0b57] | 920 | /** Prints molecule to *out.
|
---|
| 921 | * \param *out output stream
|
---|
[14de469] | 922 | */
|
---|
[e4afb4] | 923 | bool molecule::Output(ostream * const output) const
|
---|
[14de469] | 924 | {
|
---|
[e138de] | 925 | if (output == NULL) {
|
---|
[cee0b57] | 926 | return false;
|
---|
| 927 | } else {
|
---|
[0ba410] | 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;
|
---|
[30c753] | 932 | for_each(begin(),end(),boost::bind(&atom::OutputArrayIndexed,_1,output,elementLookup,AtomNo,(const char*)0));
|
---|
[cee0b57] | 933 | return true;
|
---|
[042f82] | 934 | }
|
---|
[14de469] | 935 | };
|
---|
| 936 |
|
---|
[266237] | 937 | /** Outputs contents of each atom::ListOfBonds.
|
---|
[cee0b57] | 938 | * \param *out output stream
|
---|
[14de469] | 939 | */
|
---|
[e138de] | 940 | void molecule::OutputListOfBonds() const
|
---|
[14de469] | 941 | {
|
---|
[4b5cf8] | 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 | }
|
---|
[14de469] | 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 | */
|
---|
[f01769] | 960 | int * molecule::GetFatherSonAtomicMap(const molecule * const OtherMolecule)
|
---|
[14de469] | 961 | {
|
---|
[47d041] | 962 | LOG(3, "Begin of GetFatherAtomicMap.");
|
---|
[1024cb] | 963 | int *AtomicMap = new int[getAtomCount()];
|
---|
[ea7176] | 964 | for (int i=getAtomCount();i--;)
|
---|
[042f82] | 965 | AtomicMap[i] = -1;
|
---|
| 966 | if (OtherMolecule == this) { // same molecule
|
---|
[ea7176] | 967 | for (int i=getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence
|
---|
[042f82] | 968 | AtomicMap[i] = i;
|
---|
[47d041] | 969 | LOG(4, "Map is trivial.");
|
---|
[042f82] | 970 | } else {
|
---|
[47d041] | 971 | std::stringstream output;
|
---|
| 972 | output << "Map is ";
|
---|
[f01769] | 973 | for (molecule::const_iterator iter = const_cast<const molecule &>(*this).begin();
|
---|
| 974 | iter != const_cast<const molecule &>(*this).end();
|
---|
| 975 | ++iter) {
|
---|
[910a5d] | 976 | if ((*iter)->getFather() == NULL) {
|
---|
[735b1c] | 977 | AtomicMap[(*iter)->getNr()] = -2;
|
---|
[042f82] | 978 | } else {
|
---|
[9879f6] | 979 | for (molecule::const_iterator runner = OtherMolecule->begin(); runner != OtherMolecule->end(); ++runner) {
|
---|
[042f82] | 980 | //for (int i=0;i<AtomCount;i++) { // search atom
|
---|
[1024cb] | 981 | //for (int j=0;j<OtherMolecule->getAtomCount();j++) {
|
---|
[910a5d] | 982 | //LOG(4, "Comparing father " << (*iter)->getFather() << " with the other one " << (*runner)->getFather() << ".");
|
---|
| 983 | if ((*iter)->getFather() == (*runner))
|
---|
[735b1c] | 984 | AtomicMap[(*iter)->getNr()] = (*runner)->getNr();
|
---|
[042f82] | 985 | }
|
---|
| 986 | }
|
---|
[47d041] | 987 | output << AtomicMap[(*iter)->getNr()] << "\t";
|
---|
[042f82] | 988 | }
|
---|
[47d041] | 989 | LOG(4, output.str());
|
---|
[042f82] | 990 | }
|
---|
[47d041] | 991 | LOG(3, "End of GetFatherAtomicMap.");
|
---|
[042f82] | 992 | return AtomicMap;
|
---|
[14de469] | 993 | };
|
---|
| 994 |
|
---|
[4a7776a] | 995 |
|
---|
[c68025] | 996 | void molecule::flipActiveFlag(){
|
---|
| 997 | ActiveFlag = !ActiveFlag;
|
---|
| 998 | }
|
---|
[560bbe] | 999 |
|
---|
[aeb694] | 1000 | Shape 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 |
|
---|
| 1018 | Shape molecule::getBoundingSphere(const double boundary) const
|
---|
[c67ff9] | 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();
|
---|
[8c001a] | 1027 | if (begin() != end())
|
---|
| 1028 | center *= 1./(double)size();
|
---|
[c67ff9] | 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
|
---|
[55feea1] | 1037 | radius = sqrt(radius) + boundary + 1e+6*std::numeric_limits<double>::epsilon();
|
---|
[c67ff9] | 1038 | LOG(1, "INFO: The " << size() << " atoms of the molecule are contained in a sphere at "
|
---|
| 1039 | << center << " with radius " << radius << ".");
|
---|
| 1040 |
|
---|
[f24af7] | 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.
|
---|
[c67ff9] | 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 |
|
---|
[e39e7a] | 1049 | molecule::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 |
|
---|
[cbd409] | 1065 | Vector molecule::updateMoleculeCenter() const
|
---|
| 1066 | {
|
---|
| 1067 | return (1./(double)getAtomCount())*molcenter;
|
---|
| 1068 | }
|
---|
| 1069 |
|
---|
[e39e7a] | 1070 | molecule::BoundingBoxInfo molecule::getBoundingBox() const
|
---|
| 1071 | {
|
---|
| 1072 | return **BoundingBox;
|
---|
| 1073 | }
|
---|
| 1074 |
|
---|
[cbd409] | 1075 | Vector molecule::getMoleculeCenter() const
|
---|
| 1076 | {
|
---|
| 1077 | return **MoleculeCenter;
|
---|
| 1078 | }
|
---|
| 1079 |
|
---|
[c32d21] | 1080 | void molecule::update(Observable *publisher)
|
---|
| 1081 | {
|
---|
| 1082 | ASSERT(0, "molecule::update() - did not sign on for any general updates.");
|
---|
| 1083 | }
|
---|
| 1084 |
|
---|
| 1085 | void 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
|
---|
[fb95a5] | 1098 | _lastchangedatomid = _atom->getId();
|
---|
[cbd409] | 1099 | // update entry in map and also molecule center
|
---|
[e39e7a] | 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.");
|
---|
[cbd409] | 1106 | molcenter[i] -= iter->second;
|
---|
[e39e7a] | 1107 | BoundingBoxSweepingAxis[i].left.erase(iter);
|
---|
[cbd409] | 1108 | const Vector &position = _atom->getPosition();
|
---|
[e39e7a] | 1109 | BoundingBoxSweepingAxis[i].left.insert(
|
---|
[cbd409] | 1110 | std::make_pair(_atom->getId(), position[i]) );
|
---|
| 1111 | molcenter[i] += position[i];
|
---|
[e39e7a] | 1112 | }
|
---|
| 1113 | BoundingBoxInfo newinfo = updateBoundingBox();
|
---|
[c32d21] | 1114 | OBSERVE;
|
---|
| 1115 | NOTIFY(AtomMoved);
|
---|
[cbd409] | 1116 | NOTIFY(MoleculeCenterChanged);
|
---|
[e39e7a] | 1117 | if (oldinfo != newinfo)
|
---|
| 1118 | NOTIFY(BoundingBoxChanged);
|
---|
[c32d21] | 1119 | break;
|
---|
| 1120 | }
|
---|
[9b3262b] | 1121 | case AtomObservable::ElementChanged:
|
---|
| 1122 | {
|
---|
| 1123 | // emit others about one of our atoms moved
|
---|
[fb95a5] | 1124 | _lastchangedatomid = _atom->getId();
|
---|
[9b3262b] | 1125 | OBSERVE;
|
---|
| 1126 | NOTIFY(FormulaChanged);
|
---|
[29f7c1] | 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;
|
---|
[9b3262b] | 1138 | break;
|
---|
| 1139 | }
|
---|
[29f7c1] | 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 | }
|
---|
[c32d21] | 1154 | default:
|
---|
| 1155 | ASSERT( 0, "molecule::recieveNotification() - we did not sign up for channel "
|
---|
| 1156 | +toString(notification->getChannelNo()));
|
---|
| 1157 | break;
|
---|
| 1158 | }
|
---|
| 1159 | }
|
---|
| 1160 | }
|
---|
| 1161 |
|
---|
| 1162 | void molecule::subjectKilled(Observable *publisher)
|
---|
| 1163 | {
|
---|
| 1164 | // do nothing, atom does it all
|
---|
| 1165 | }
|
---|
| 1166 |
|
---|
[b71881] | 1167 | void molecule::select()
|
---|
| 1168 | {
|
---|
| 1169 | OBSERVE;
|
---|
| 1170 | selected = true;
|
---|
| 1171 | NOTIFY(SelectionChanged);
|
---|
| 1172 | }
|
---|
| 1173 |
|
---|
| 1174 | void molecule::unselect()
|
---|
| 1175 | {
|
---|
| 1176 | OBSERVE;
|
---|
| 1177 | selected = false;
|
---|
| 1178 | NOTIFY(SelectionChanged);
|
---|
| 1179 | }
|
---|
| 1180 |
|
---|
[c0f2fc] | 1181 | void molecule::associateAtomWithMolecule(atom *_atom)
|
---|
| 1182 | {
|
---|
| 1183 | _atom->signOn(this, AtomObservable::PositionChanged);
|
---|
| 1184 | _atom->signOn(this, AtomObservable::ElementChanged);
|
---|
[29f7c1] | 1185 | _atom->signOn(this, AtomObservable::BondsAdded);
|
---|
| 1186 | _atom->signOn(this, AtomObservable::BondsRemoved);
|
---|
[c0f2fc] | 1187 | insert(_atom);
|
---|
[29f7c1] | 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 | {
|
---|
[d316ad] | 1200 | OBSERVE;
|
---|
[29f7c1] | 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;
|
---|
[d316ad] | 1210 | NOTIFY(FormulaChanged);
|
---|
[29f7c1] | 1211 | }
|
---|
[c0f2fc] | 1212 | }
|
---|
| 1213 |
|
---|
| 1214 | void molecule::disassociateAtomWithMolecule(atom *_atom)
|
---|
| 1215 | {
|
---|
| 1216 | _atom->signOff(this, AtomObservable::PositionChanged);
|
---|
| 1217 | _atom->signOff(this, AtomObservable::ElementChanged);
|
---|
[29f7c1] | 1218 | _atom->signOff(this, AtomObservable::BondsAdded);
|
---|
| 1219 | _atom->signOff(this, AtomObservable::BondsRemoved);
|
---|
[c0f2fc] | 1220 | erase(_atom);
|
---|
[29f7c1] | 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 | {
|
---|
[d316ad] | 1230 | OBSERVE;
|
---|
[29f7c1] | 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);
|
---|
[d316ad] | 1237 | NOTIFY(FormulaChanged);
|
---|
[29f7c1] | 1238 | }
|
---|
[c0f2fc] | 1239 | }
|
---|
| 1240 |
|
---|
[560bbe] | 1241 | // construct idpool
|
---|
| 1242 | CONSTRUCT_IDPOOL(atomId_t, continuousId)
|
---|
[c67ff9] | 1243 |
|
---|