| 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 atom.cpp
 | 
|---|
| 25 |  *
 | 
|---|
| 26 |  * Function implementations for the class atom.
 | 
|---|
| 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 "atom.hpp"
 | 
|---|
| 38 | #include "AtomObserver.hpp"
 | 
|---|
| 39 | #include "Bond/bond.hpp"
 | 
|---|
| 40 | #include "CodePatterns/Log.hpp"
 | 
|---|
| 41 | #include "config.hpp"
 | 
|---|
| 42 | #include "Element/element.hpp"
 | 
|---|
| 43 | #include "LinearAlgebra/Vector.hpp"
 | 
|---|
| 44 | #include "World.hpp"
 | 
|---|
| 45 | #include "WorldTime.hpp"
 | 
|---|
| 46 | #include "molecule.hpp"
 | 
|---|
| 47 | #include "Shapes/Shape.hpp"
 | 
|---|
| 48 | 
 | 
|---|
| 49 | #include <iomanip>
 | 
|---|
| 50 | #include <iostream>
 | 
|---|
| 51 | 
 | 
|---|
| 52 | /************************************* Functions for class atom *************************************/
 | 
|---|
| 53 | 
 | 
|---|
| 54 | 
 | 
|---|
| 55 | atom::atom() :
 | 
|---|
| 56 |   father(this),
 | 
|---|
| 57 |   sort(&Nr),
 | 
|---|
| 58 |   mol(0),
 | 
|---|
| 59 |   selected(false)
 | 
|---|
| 60 | {
 | 
|---|
| 61 |   // note AtomObserver about inserted atom
 | 
|---|
| 62 |   AtomObserver::getInstance().Inserted(this);
 | 
|---|
| 63 | }
 | 
|---|
| 64 | 
 | 
|---|
| 65 | atom::atom(atom *pointer) :
 | 
|---|
| 66 |     ParticleInfo(*pointer),
 | 
|---|
| 67 |     AtomInfo(*pointer),
 | 
|---|
| 68 |     father(pointer),
 | 
|---|
| 69 |     sort(&Nr),
 | 
|---|
| 70 |     mol(0),
 | 
|---|
| 71 |     selected(false)
 | 
|---|
| 72 | {
 | 
|---|
| 73 |   // sign on to father atom to be notified when it is removed
 | 
|---|
| 74 |   father->signOn(this);
 | 
|---|
| 75 | 
 | 
|---|
| 76 |   // note AtomObserver about inserted atom
 | 
|---|
| 77 |   AtomObserver::getInstance().Inserted(this);
 | 
|---|
| 78 | };
 | 
|---|
| 79 | 
 | 
|---|
| 80 | atom *atom::clone(){
 | 
|---|
| 81 |   atom *res = new atom(this);
 | 
|---|
| 82 |   World::getInstance().registerAtom(res);
 | 
|---|
| 83 |   return res;
 | 
|---|
| 84 | }
 | 
|---|
| 85 | 
 | 
|---|
| 86 | 
 | 
|---|
| 87 | /** Destructor of class atom.
 | 
|---|
| 88 |  */
 | 
|---|
| 89 | atom::~atom()
 | 
|---|
| 90 | {
 | 
|---|
| 91 |   // sign off from possible father
 | 
|---|
| 92 |   if ((father != this) && (father != NULL))
 | 
|---|
| 93 |     father->signOff(this);
 | 
|---|
| 94 | 
 | 
|---|
| 95 |   removeFromMolecule();
 | 
|---|
| 96 |   // note AtomObserver about removed atom
 | 
|---|
| 97 |   AtomObserver::getInstance().Removed(this);
 | 
|---|
| 98 | }
 | 
|---|
| 99 | 
 | 
|---|
| 100 | 
 | 
|---|
| 101 | void atom::UpdateStep(const unsigned int _step)
 | 
|---|
| 102 | {
 | 
|---|
| 103 |   LOG(4,"atom::UpdateStep() called.");
 | 
|---|
| 104 |   // append to position, velocity and force vector
 | 
|---|
| 105 |   AtomInfo::AppendTrajectoryStep(WorldTime::getTime());
 | 
|---|
| 106 |   // append to ListOfBonds vector
 | 
|---|
| 107 |   BondedParticleInfo::AppendTrajectoryStep(WorldTime::getTime());
 | 
|---|
| 108 | }
 | 
|---|
| 109 | 
 | 
|---|
| 110 | void atom::removeSteps(const unsigned int _firststep, const unsigned int _laststep)
 | 
|---|
| 111 | {
 | 
|---|
| 112 |   LOG(4,"atom::removeSteps() called.");
 | 
|---|
| 113 |   // append to position, velocity and force vector
 | 
|---|
| 114 |   AtomInfo::removeTrajectorySteps(_firststep, _laststep);
 | 
|---|
| 115 |   // append to ListOfBonds vector
 | 
|---|
| 116 |   BondedParticleInfo::removeTrajectorySteps(_firststep, _laststep);
 | 
|---|
| 117 | }
 | 
|---|
| 118 | 
 | 
|---|
| 119 | atom *atom::GetTrueFather()
 | 
|---|
| 120 | {
 | 
|---|
| 121 |   const atom *father = const_cast<const atom *>(this)->GetTrueFather();
 | 
|---|
| 122 |   return const_cast<atom *>(father);
 | 
|---|
| 123 | }
 | 
|---|
| 124 | 
 | 
|---|
| 125 | const atom *atom::GetTrueFather() const
 | 
|---|
| 126 | {
 | 
|---|
| 127 |   if(father == this){ // top most father is the one that points on itself
 | 
|---|
| 128 |     return this;
 | 
|---|
| 129 |   }
 | 
|---|
| 130 |   else if(!father) {
 | 
|---|
| 131 |     return 0;
 | 
|---|
| 132 |   }
 | 
|---|
| 133 |   else {
 | 
|---|
| 134 |     return father->GetTrueFather();
 | 
|---|
| 135 |   }
 | 
|---|
| 136 | }
 | 
|---|
| 137 | 
 | 
|---|
| 138 | void atom::setFather(atom * const _father)
 | 
|---|
| 139 | {
 | 
|---|
| 140 |   // sign off from old father
 | 
|---|
| 141 |   if ((father != this) && (father != NULL))
 | 
|---|
| 142 |     father->signOff(this);
 | 
|---|
| 143 | 
 | 
|---|
| 144 |   father = _father;
 | 
|---|
| 145 |   father->signOn(this);
 | 
|---|
| 146 | }
 | 
|---|
| 147 | 
 | 
|---|
| 148 | /** Sets father to itself or its father in case of copying a molecule.
 | 
|---|
| 149 |  */
 | 
|---|
| 150 | void atom::CorrectFather()
 | 
|---|
| 151 | {
 | 
|---|
| 152 |   if (father->father != father)   // same atom in copy's father points to itself
 | 
|---|
| 153 | //    father = this;  // set father to itself (copy of a whole molecule)
 | 
|---|
| 154 | //  else
 | 
|---|
| 155 |    father = father->father;  // set father to original's father
 | 
|---|
| 156 | 
 | 
|---|
| 157 | };
 | 
|---|
| 158 | 
 | 
|---|
| 159 | void atom::EqualsFather ( const atom *ptr, const atom **res ) const
 | 
|---|
| 160 | {
 | 
|---|
| 161 |   if ( ptr == father )
 | 
|---|
| 162 |     *res = this;
 | 
|---|
| 163 | };
 | 
|---|
| 164 | 
 | 
|---|
| 165 | bool atom::isFather(const atom *ptr){
 | 
|---|
| 166 |   return ptr==father;
 | 
|---|
| 167 | }
 | 
|---|
| 168 | 
 | 
|---|
| 169 | bool atom::OutputIndexed(ofstream * const out, const int ElementNo, const int AtomNo, const char *comment) const
 | 
|---|
| 170 | {
 | 
|---|
| 171 |   if (out != NULL) {
 | 
|---|
| 172 |     *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t"  << fixed << setprecision(9) << showpoint;
 | 
|---|
| 173 |     *out << at(0) << "\t" << at(1) << "\t" << at(2);
 | 
|---|
| 174 |     *out << "\t" << (int)(getFixedIon());
 | 
|---|
| 175 |     if (getAtomicVelocity().Norm() > MYEPSILON)
 | 
|---|
| 176 |       *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t";
 | 
|---|
| 177 |     if (comment != NULL)
 | 
|---|
| 178 |       *out << " # " << comment << endl;
 | 
|---|
| 179 |     else
 | 
|---|
| 180 |       *out << " # molecule nr " << getNr() << endl;
 | 
|---|
| 181 |     return true;
 | 
|---|
| 182 |   } else
 | 
|---|
| 183 |     return false;
 | 
|---|
| 184 | };
 | 
|---|
| 185 | 
 | 
|---|
| 186 | bool atom::OutputArrayIndexed(ostream * const out,const enumeration<const element*> &elementLookup, int *AtomNo, const char *comment) const
 | 
|---|
| 187 | {
 | 
|---|
| 188 |   AtomNo[getType()->getAtomicNumber()]++;  // increment number
 | 
|---|
| 189 |   if (out != NULL) {
 | 
|---|
| 190 |     const element *elemental = getType();
 | 
|---|
| 191 |     ASSERT(elementLookup.there.find(elemental)!=elementLookup.there.end(),"Type of this atom was not in the formula upon enumeration");
 | 
|---|
| 192 |     *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[elemental->getAtomicNumber()] << "\t"  << fixed << setprecision(9) << showpoint;
 | 
|---|
| 193 |     *out << at(0) << "\t" << at(1) << "\t" << at(2);
 | 
|---|
| 194 |     *out << "\t" << getFixedIon();
 | 
|---|
| 195 |     if (getAtomicVelocity().Norm() > MYEPSILON)
 | 
|---|
| 196 |       *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t";
 | 
|---|
| 197 |     if (comment != NULL)
 | 
|---|
| 198 |       *out << " # " << comment << endl;
 | 
|---|
| 199 |     else
 | 
|---|
| 200 |       *out << " # molecule nr " << getNr() << endl;
 | 
|---|
| 201 |     return true;
 | 
|---|
| 202 |   } else
 | 
|---|
| 203 |     return false;
 | 
|---|
| 204 | };
 | 
|---|
| 205 | 
 | 
|---|
| 206 | bool atom::Compare(const atom &ptr) const
 | 
|---|
| 207 | {
 | 
|---|
| 208 |   if (getNr() < ptr.getNr())
 | 
|---|
| 209 |     return true;
 | 
|---|
| 210 |   else
 | 
|---|
| 211 |     return false;
 | 
|---|
| 212 | };
 | 
|---|
| 213 | 
 | 
|---|
| 214 | double atom::DistanceSquaredToVector(const Vector &origin) const
 | 
|---|
| 215 | {
 | 
|---|
| 216 |   return DistanceSquared(origin);
 | 
|---|
| 217 | };
 | 
|---|
| 218 | 
 | 
|---|
| 219 | double atom::DistanceToVector(const Vector &origin) const
 | 
|---|
| 220 | {
 | 
|---|
| 221 |   return distance(origin);
 | 
|---|
| 222 | };
 | 
|---|
| 223 | 
 | 
|---|
| 224 | void atom::InitComponentNr()
 | 
|---|
| 225 | {
 | 
|---|
| 226 |   if (ComponentNr != NULL)
 | 
|---|
| 227 |     delete[](ComponentNr);
 | 
|---|
| 228 |   const BondList& ListOfBonds = getListOfBonds();
 | 
|---|
| 229 |   ComponentNr = new int[ListOfBonds.size()+1];
 | 
|---|
| 230 |   for (int i=ListOfBonds.size()+1;i--;)
 | 
|---|
| 231 |     ComponentNr[i] = -1;
 | 
|---|
| 232 | };
 | 
|---|
| 233 | 
 | 
|---|
| 234 | void atom::resetGraphNr(){
 | 
|---|
| 235 |   GraphNr=-1;
 | 
|---|
| 236 | }
 | 
|---|
| 237 | 
 | 
|---|
| 238 | std::ostream & atom::operator << (std::ostream &ost) const
 | 
|---|
| 239 | {
 | 
|---|
| 240 |   ParticleInfo::operator<<(ost);
 | 
|---|
| 241 |   ost << "," << getPosition();
 | 
|---|
| 242 |   return ost;
 | 
|---|
| 243 | }
 | 
|---|
| 244 | 
 | 
|---|
| 245 | std::ostream & operator << (std::ostream &ost, const atom &a)
 | 
|---|
| 246 | {
 | 
|---|
| 247 |   a.ParticleInfo::operator<<(ost);
 | 
|---|
| 248 |   ost << "," << a.getPosition();
 | 
|---|
| 249 |   return ost;
 | 
|---|
| 250 | }
 | 
|---|
| 251 | 
 | 
|---|
| 252 | bool operator < (atom &a, atom &b)
 | 
|---|
| 253 | {
 | 
|---|
| 254 |   return a.Compare(b);
 | 
|---|
| 255 | };
 | 
|---|
| 256 | 
 | 
|---|
| 257 | World *atom::getWorld(){
 | 
|---|
| 258 |   return world;
 | 
|---|
| 259 | }
 | 
|---|
| 260 | 
 | 
|---|
| 261 | void atom::setWorld(World* _world){
 | 
|---|
| 262 |   world = _world;
 | 
|---|
| 263 | }
 | 
|---|
| 264 | 
 | 
|---|
| 265 | bool atom::changeId(atomId_t newId){
 | 
|---|
| 266 |   // first we move ourselves in the world
 | 
|---|
| 267 |   // the world lets us know if that succeeded
 | 
|---|
| 268 |   atomId_t oldid = id;
 | 
|---|
| 269 |   if(world->changeAtomId(id,newId,this)){
 | 
|---|
| 270 |     OBSERVE;
 | 
|---|
| 271 |     id = newId;
 | 
|---|
| 272 |     if (mol != NULL)
 | 
|---|
| 273 |       mol->changeAtomId(oldid, newId);
 | 
|---|
| 274 |     NOTIFY(IndexChanged);
 | 
|---|
| 275 |     return true;
 | 
|---|
| 276 |   }
 | 
|---|
| 277 |   else{
 | 
|---|
| 278 |     return false;
 | 
|---|
| 279 |   }
 | 
|---|
| 280 | }
 | 
|---|
| 281 | 
 | 
|---|
| 282 | void atom::setId(atomId_t _id) {
 | 
|---|
| 283 |   id=_id;
 | 
|---|
| 284 | }
 | 
|---|
| 285 | 
 | 
|---|
| 286 | atomId_t atom::getId() const {
 | 
|---|
| 287 |   return id;
 | 
|---|
| 288 | }
 | 
|---|
| 289 | 
 | 
|---|
| 290 | void atom::setMolecule(molecule *_mol){
 | 
|---|
| 291 |   {
 | 
|---|
| 292 |     // we must inform about molecule changed before the molecule itself tells about it
 | 
|---|
| 293 |     OBSERVE;
 | 
|---|
| 294 |     NOTIFY(MoleculeChanged);
 | 
|---|
| 295 |     // take this atom from the old molecule
 | 
|---|
| 296 |     removeFromMolecule();
 | 
|---|
| 297 |     mol = _mol;
 | 
|---|
| 298 |   }
 | 
|---|
| 299 |   if ((mol) && (!mol->containsAtom(this))) {
 | 
|---|
| 300 |     mol->associateAtomWithMolecule(this);
 | 
|---|
| 301 |   }
 | 
|---|
| 302 | }
 | 
|---|
| 303 | 
 | 
|---|
| 304 | const molecule* atom::getMolecule() const {
 | 
|---|
| 305 |   return mol;
 | 
|---|
| 306 | }
 | 
|---|
| 307 | 
 | 
|---|
| 308 | void atom::removeFromMolecule(){
 | 
|---|
| 309 |   if(mol){
 | 
|---|
| 310 |     OBSERVE;
 | 
|---|
| 311 |     NOTIFY(MoleculeChanged);
 | 
|---|
| 312 |     if(mol->containsAtom(this)){
 | 
|---|
| 313 |       mol->disassociateAtomWithMolecule(this);
 | 
|---|
| 314 |     }
 | 
|---|
| 315 |     mol=0;
 | 
|---|
| 316 |   }
 | 
|---|
| 317 | }
 | 
|---|
| 318 | 
 | 
|---|
| 319 | bool atom::changeNr(const int newNr)
 | 
|---|
| 320 | {
 | 
|---|
| 321 |   if ((mol) && (mol->changeAtomNr(getNr(),newNr,this))) {
 | 
|---|
| 322 |     return true;
 | 
|---|
| 323 |   } else{
 | 
|---|
| 324 |     return false;
 | 
|---|
| 325 |   }
 | 
|---|
| 326 | }
 | 
|---|
| 327 | 
 | 
|---|
| 328 | int atom::getNr() const{
 | 
|---|
| 329 |   return ParticleInfo::getNr();
 | 
|---|
| 330 | }
 | 
|---|
| 331 | 
 | 
|---|
| 332 | atom* NewAtom(atomId_t _id){
 | 
|---|
| 333 |   atom * res = new atom();
 | 
|---|
| 334 |   res->setId(_id);
 | 
|---|
| 335 |   return res;
 | 
|---|
| 336 | }
 | 
|---|
| 337 | 
 | 
|---|
| 338 | void DeleteAtom(atom* atom){
 | 
|---|
| 339 |   delete atom;
 | 
|---|
| 340 | }
 | 
|---|
| 341 | 
 | 
|---|
| 342 | bool compareAtomElements(atom* atom1,atom* atom2){
 | 
|---|
| 343 |   return atom1->getType()->getAtomicNumber() < atom2->getType()->getAtomicNumber();
 | 
|---|
| 344 | }
 | 
|---|
| 345 | /*
 | 
|---|
| 346 | void atom::update(Observable *publisher)
 | 
|---|
| 347 | {}
 | 
|---|
| 348 | 
 | 
|---|
| 349 | void atom::recieveNotification(Observable *publisher, Notification_ptr notification)
 | 
|---|
| 350 | {
 | 
|---|
| 351 |   ASSERT(0, "atom::recieveNotification() - we are not signed on to any notifications.");
 | 
|---|
| 352 | }
 | 
|---|
| 353 | */
 | 
|---|
| 354 | void atom::subjectKilled(Observable *publisher)
 | 
|---|
| 355 | {
 | 
|---|
| 356 |   // as publisher has been half-deallocated (Observable is one of the base classes, hence
 | 
|---|
| 357 |   // becomes destroyed latest), we cannot senibly cast it anymore.
 | 
|---|
| 358 |   // Hence, we simply have to check here whether it is NOT one of the other instances
 | 
|---|
| 359 |   // we are signed on to.
 | 
|---|
| 360 |   father = this;
 | 
|---|
| 361 |   // no need to sign off
 | 
|---|
| 362 | }
 | 
|---|
| 363 | 
 | 
|---|
| 364 | void atom::select()
 | 
|---|
| 365 | {
 | 
|---|
| 366 |   OBSERVE;
 | 
|---|
| 367 |   selected = true;
 | 
|---|
| 368 |   NOTIFY(SelectionChanged);
 | 
|---|
| 369 | }
 | 
|---|
| 370 | 
 | 
|---|
| 371 | void atom::unselect()
 | 
|---|
| 372 | {
 | 
|---|
| 373 |   OBSERVE;
 | 
|---|
| 374 |   selected = false;
 | 
|---|
| 375 |   NOTIFY(SelectionChanged);
 | 
|---|
| 376 | }
 | 
|---|