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