| [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 atom.cpp
 | 
|---|
| [1907a7] | 25 |  *
 | 
|---|
| [14de469] | 26 |  * Function implementations for the class atom.
 | 
|---|
| [1907a7] | 27 |  *
 | 
|---|
| [14de469] | 28 |  */
 | 
|---|
 | 29 | 
 | 
|---|
| [bf3817] | 30 | // include config.h
 | 
|---|
 | 31 | #ifdef HAVE_CONFIG_H
 | 
|---|
 | 32 | #include <config.h>
 | 
|---|
 | 33 | #endif
 | 
|---|
 | 34 | 
 | 
|---|
| [ad011c] | 35 | #include "CodePatterns/MemDebug.hpp"
 | 
|---|
| [112b09] | 36 | 
 | 
|---|
| [357fba] | 37 | #include "atom.hpp"
 | 
|---|
| [708277] | 38 | #include "AtomObserver.hpp"
 | 
|---|
| [129204] | 39 | #include "Bond/bond.hpp"
 | 
|---|
| [e2373df] | 40 | #include "CodePatterns/Log.hpp"
 | 
|---|
| [4a7776a] | 41 | #include "config.hpp"
 | 
|---|
| [3bdb6d] | 42 | #include "Element/element.hpp"
 | 
|---|
| [57f243] | 43 | #include "LinearAlgebra/Vector.hpp"
 | 
|---|
| [d346b6] | 44 | #include "World.hpp"
 | 
|---|
| [11f0fa] | 45 | #include "WorldTime.hpp"
 | 
|---|
| [6cfa36] | 46 | #include "molecule.hpp"
 | 
|---|
| [c550dd] | 47 | #include "Shapes/Shape.hpp"
 | 
|---|
| [a0064e] | 48 | 
 | 
|---|
| [36166d] | 49 | #include <iomanip>
 | 
|---|
| [0ba410] | 50 | #include <iostream>
 | 
|---|
| [36166d] | 51 | 
 | 
|---|
| [14de469] | 52 | /************************************* Functions for class atom *************************************/
 | 
|---|
 | 53 | 
 | 
|---|
| [70ff32] | 54 | 
 | 
|---|
| [d05088] | 55 | atom::atom() :
 | 
|---|
| [97b825] | 56 |   father(this),
 | 
|---|
| [5309ba] | 57 |   sort(&Nr),
 | 
|---|
| [97b825] | 58 |   mol(0)
 | 
|---|
| [708277] | 59 | {
 | 
|---|
| [d05088] | 60 |   // note AtomObserver about inserted atom
 | 
|---|
| [708277] | 61 |   AtomObserver::getInstance().AtomInserted(this);
 | 
|---|
 | 62 | }
 | 
|---|
| [14de469] | 63 | 
 | 
|---|
| [ea0c8b] | 64 | atom::atom(atom *pointer) :
 | 
|---|
 | 65 |     ParticleInfo(*pointer),
 | 
|---|
| [c742fe] | 66 |     AtomInfo(*pointer),
 | 
|---|
| [97b825] | 67 |     father(pointer),
 | 
|---|
| [5309ba] | 68 |     sort(&Nr),
 | 
|---|
| [9df680] | 69 |     mol(0)
 | 
|---|
| [2319ed] | 70 | {
 | 
|---|
| [d05088] | 71 |   // sign on to father atom to be notified when it is removed
 | 
|---|
 | 72 |   father->signOn(this);
 | 
|---|
 | 73 | 
 | 
|---|
 | 74 |   // note AtomObserver about inserted atom
 | 
|---|
| [708277] | 75 |   AtomObserver::getInstance().AtomInserted(this);
 | 
|---|
| [b453f9] | 76 | };
 | 
|---|
| [2319ed] | 77 | 
 | 
|---|
| [46d958] | 78 | atom *atom::clone(){
 | 
|---|
| [68f03d] | 79 |   atom *res = new atom(this);
 | 
|---|
| [23b547] | 80 |   World::getInstance().registerAtom(res);
 | 
|---|
| [46d958] | 81 |   return res;
 | 
|---|
 | 82 | }
 | 
|---|
 | 83 | 
 | 
|---|
| [2319ed] | 84 | 
 | 
|---|
| [14de469] | 85 | /** Destructor of class atom.
 | 
|---|
 | 86 |  */
 | 
|---|
| [1907a7] | 87 | atom::~atom()
 | 
|---|
| [14de469] | 88 | {
 | 
|---|
| [d05088] | 89 |   // sign off from possible father
 | 
|---|
 | 90 |   if ((father != this) && (father != NULL))
 | 
|---|
 | 91 |     father->signOff(this);
 | 
|---|
 | 92 | 
 | 
|---|
| [6cfa36] | 93 |   removeFromMolecule();
 | 
|---|
| [d05088] | 94 |   // note AtomObserver about removed atom
 | 
|---|
| [708277] | 95 |   AtomObserver::getInstance().AtomRemoved(this);
 | 
|---|
 | 96 | }
 | 
|---|
| [e2373df] | 97 | 
 | 
|---|
 | 98 | 
 | 
|---|
| [8cc22f] | 99 | void atom::UpdateStep(const unsigned int _step)
 | 
|---|
| [e2373df] | 100 | {
 | 
|---|
| [8cc22f] | 101 |   LOG(4,"atom::UpdateStep() called.");
 | 
|---|
| [e2373df] | 102 |   // append to position, velocity and force vector
 | 
|---|
| [8cc22f] | 103 |   AtomInfo::AppendTrajectoryStep(WorldTime::getTime()+1);
 | 
|---|
| [1e6249] | 104 |   // append to ListOfBonds vector
 | 
|---|
| [8cc22f] | 105 |   BondedParticleInfo::AppendTrajectoryStep(WorldTime::getTime()+1);
 | 
|---|
| [e2373df] | 106 | }
 | 
|---|
 | 107 | 
 | 
|---|
| [8cc22f] | 108 | void atom::removeStep(const unsigned int _step)
 | 
|---|
| [7e51e1] | 109 | {
 | 
|---|
| [8cc22f] | 110 |   LOG(4,"atom::removeStep() called.");
 | 
|---|
| [7e51e1] | 111 |   // append to position, velocity and force vector
 | 
|---|
| [8cc22f] | 112 |   AtomInfo::removeTrajectoryStep(_step);
 | 
|---|
| [7e51e1] | 113 |   // append to ListOfBonds vector
 | 
|---|
| [8cc22f] | 114 |   BondedParticleInfo::removeTrajectoryStep(_step);
 | 
|---|
| [7e51e1] | 115 | }
 | 
|---|
 | 116 | 
 | 
|---|
| [59fff1] | 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
 | 
|---|
| [14de469] | 124 | {
 | 
|---|
| [215df0] | 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 |   }
 | 
|---|
| [910a5d] | 134 | }
 | 
|---|
 | 135 | 
 | 
|---|
 | 136 | void atom::setFather(atom * const _father)
 | 
|---|
 | 137 | {
 | 
|---|
| [d05088] | 138 |   // sign off from old father
 | 
|---|
 | 139 |   if ((father != this) && (father != NULL))
 | 
|---|
 | 140 |     father->signOff(this);
 | 
|---|
 | 141 | 
 | 
|---|
| [910a5d] | 142 |   father = _father;
 | 
|---|
| [d05088] | 143 |   father->signOn(this);
 | 
|---|
| [910a5d] | 144 | }
 | 
|---|
| [14de469] | 145 | 
 | 
|---|
| [e65246] | 146 | /** Sets father to itself or its father in case of copying a molecule.
 | 
|---|
 | 147 |  */
 | 
|---|
 | 148 | void atom::CorrectFather()
 | 
|---|
 | 149 | {
 | 
|---|
| [2e352f] | 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
 | 
|---|
| [e65246] | 153 |    father = father->father;  // set father to original's father
 | 
|---|
 | 154 | 
 | 
|---|
 | 155 | };
 | 
|---|
 | 156 | 
 | 
|---|
| [b453f9] | 157 | void atom::EqualsFather ( const atom *ptr, const atom **res ) const
 | 
|---|
| [e65246] | 158 | {
 | 
|---|
 | 159 |   if ( ptr == father )
 | 
|---|
 | 160 |     *res = this;
 | 
|---|
 | 161 | };
 | 
|---|
 | 162 | 
 | 
|---|
| [00abfc] | 163 | bool atom::isFather(const atom *ptr){
 | 
|---|
 | 164 |   return ptr==father;
 | 
|---|
 | 165 | }
 | 
|---|
 | 166 | 
 | 
|---|
| [e138de] | 167 | bool atom::OutputIndexed(ofstream * const out, const int ElementNo, const int AtomNo, const char *comment) const
 | 
|---|
| [14de469] | 168 | {
 | 
|---|
 | 169 |   if (out != NULL) {
 | 
|---|
 | 170 |     *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t"  << fixed << setprecision(9) << showpoint;
 | 
|---|
| [d74077] | 171 |     *out << at(0) << "\t" << at(1) << "\t" << at(2);
 | 
|---|
| [6625c3] | 172 |     *out << "\t" << (int)(getFixedIon());
 | 
|---|
| [bce72c] | 173 |     if (getAtomicVelocity().Norm() > MYEPSILON)
 | 
|---|
 | 174 |       *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t";
 | 
|---|
| [437922] | 175 |     if (comment != NULL)
 | 
|---|
 | 176 |       *out << " # " << comment << endl;
 | 
|---|
| [e9f8f9] | 177 |     else
 | 
|---|
| [735b1c] | 178 |       *out << " # molecule nr " << getNr() << endl;
 | 
|---|
| [e9f8f9] | 179 |     return true;
 | 
|---|
 | 180 |   } else
 | 
|---|
 | 181 |     return false;
 | 
|---|
 | 182 | };
 | 
|---|
| [b453f9] | 183 | 
 | 
|---|
| [0ba410] | 184 | bool atom::OutputArrayIndexed(ostream * const out,const enumeration<const element*> &elementLookup, int *AtomNo, const char *comment) const
 | 
|---|
| [e9f8f9] | 185 | {
 | 
|---|
| [83f176] | 186 |   AtomNo[getType()->getAtomicNumber()]++;  // increment number
 | 
|---|
| [e9f8f9] | 187 |   if (out != NULL) {
 | 
|---|
| [8f4df1] | 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");
 | 
|---|
| [83f176] | 190 |     *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[elemental->getAtomicNumber()] << "\t"  << fixed << setprecision(9) << showpoint;
 | 
|---|
| [d74077] | 191 |     *out << at(0) << "\t" << at(1) << "\t" << at(2);
 | 
|---|
| [6625c3] | 192 |     *out << "\t" << getFixedIon();
 | 
|---|
| [bce72c] | 193 |     if (getAtomicVelocity().Norm() > MYEPSILON)
 | 
|---|
 | 194 |       *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t";
 | 
|---|
| [e9f8f9] | 195 |     if (comment != NULL)
 | 
|---|
 | 196 |       *out << " # " << comment << endl;
 | 
|---|
| [437922] | 197 |     else
 | 
|---|
| [735b1c] | 198 |       *out << " # molecule nr " << getNr() << endl;
 | 
|---|
| [14de469] | 199 |     return true;
 | 
|---|
 | 200 |   } else
 | 
|---|
 | 201 |     return false;
 | 
|---|
 | 202 | };
 | 
|---|
 | 203 | 
 | 
|---|
| [b453f9] | 204 | bool atom::Compare(const atom &ptr) const
 | 
|---|
| [4455f4] | 205 | {
 | 
|---|
| [735b1c] | 206 |   if (getNr() < ptr.getNr())
 | 
|---|
| [4455f4] | 207 |     return true;
 | 
|---|
 | 208 |   else
 | 
|---|
 | 209 |     return false;
 | 
|---|
 | 210 | };
 | 
|---|
 | 211 | 
 | 
|---|
| [b453f9] | 212 | double atom::DistanceSquaredToVector(const Vector &origin) const
 | 
|---|
| [4455f4] | 213 | {
 | 
|---|
| [d74077] | 214 |   return DistanceSquared(origin);
 | 
|---|
| [4455f4] | 215 | };
 | 
|---|
 | 216 | 
 | 
|---|
| [b453f9] | 217 | double atom::DistanceToVector(const Vector &origin) const
 | 
|---|
| [4455f4] | 218 | {
 | 
|---|
| [d74077] | 219 |   return distance(origin);
 | 
|---|
| [4455f4] | 220 | };
 | 
|---|
 | 221 | 
 | 
|---|
 | 222 | void atom::InitComponentNr()
 | 
|---|
 | 223 | {
 | 
|---|
 | 224 |   if (ComponentNr != NULL)
 | 
|---|
| [920c70] | 225 |     delete[](ComponentNr);
 | 
|---|
| [9d83b6] | 226 |   const BondList& ListOfBonds = getListOfBonds();
 | 
|---|
| [920c70] | 227 |   ComponentNr = new int[ListOfBonds.size()+1];
 | 
|---|
| [4455f4] | 228 |   for (int i=ListOfBonds.size()+1;i--;)
 | 
|---|
 | 229 |     ComponentNr[i] = -1;
 | 
|---|
| [14b65e] | 230 | };
 | 
|---|
 | 231 | 
 | 
|---|
 | 232 | void atom::resetGraphNr(){
 | 
|---|
 | 233 |   GraphNr=-1;
 | 
|---|
 | 234 | }
 | 
|---|
| [4455f4] | 235 | 
 | 
|---|
| [d74077] | 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 | }
 | 
|---|
| [4455f4] | 249 | 
 | 
|---|
 | 250 | bool operator < (atom &a, atom &b)
 | 
|---|
 | 251 | {
 | 
|---|
 | 252 |   return a.Compare(b);
 | 
|---|
 | 253 | };
 | 
|---|
 | 254 | 
 | 
|---|
| [46d958] | 255 | World *atom::getWorld(){
 | 
|---|
 | 256 |   return world;
 | 
|---|
 | 257 | }
 | 
|---|
 | 258 | 
 | 
|---|
 | 259 | void atom::setWorld(World* _world){
 | 
|---|
 | 260 |   world = _world;
 | 
|---|
 | 261 | }
 | 
|---|
 | 262 | 
 | 
|---|
| [88d586] | 263 | bool atom::changeId(atomId_t newId){
 | 
|---|
 | 264 |   // first we move ourselves in the world
 | 
|---|
 | 265 |   // the world lets us know if that succeeded
 | 
|---|
| [4f7f0bf] | 266 |   if(world->changeAtomId(id,newId,this)){
 | 
|---|
 | 267 |     OBSERVE;
 | 
|---|
 | 268 |     id = newId;
 | 
|---|
 | 269 |     NOTIFY(IndexChanged);
 | 
|---|
| [88d586] | 270 |     return true;
 | 
|---|
 | 271 |   }
 | 
|---|
 | 272 |   else{
 | 
|---|
 | 273 |     return false;
 | 
|---|
 | 274 |   }
 | 
|---|
 | 275 | }
 | 
|---|
 | 276 | 
 | 
|---|
 | 277 | void atom::setId(atomId_t _id) {
 | 
|---|
| [46d958] | 278 |   id=_id;
 | 
|---|
 | 279 | }
 | 
|---|
 | 280 | 
 | 
|---|
| [ad2b411] | 281 | atomId_t atom::getId() const {
 | 
|---|
| [46d958] | 282 |   return id;
 | 
|---|
 | 283 | }
 | 
|---|
 | 284 | 
 | 
|---|
| [6cfa36] | 285 | void atom::setMolecule(molecule *_mol){
 | 
|---|
 | 286 |   // take this atom from the old molecule
 | 
|---|
 | 287 |   removeFromMolecule();
 | 
|---|
| [3867a7] | 288 |   mol = _mol;
 | 
|---|
| [c32d21] | 289 |   if ((mol) && (!mol->containsAtom(this))) {
 | 
|---|
 | 290 |     signOn(mol, AtomObservable::PositionChanged);
 | 
|---|
 | 291 |     mol->insert(this);
 | 
|---|
 | 292 |   }
 | 
|---|
| [6cfa36] | 293 | }
 | 
|---|
 | 294 | 
 | 
|---|
| [0d9546] | 295 | void atom::unsetMolecule()
 | 
|---|
 | 296 | {
 | 
|---|
 | 297 |   // take this atom from the old molecule
 | 
|---|
 | 298 |   ASSERT(!mol->containsAtom(this),
 | 
|---|
 | 299 |       "atom::unsetMolecule() - old molecule "+toString(mol)+" still contains us!");
 | 
|---|
| [c32d21] | 300 |   signOff(mol, AtomObservable::PositionChanged);
 | 
|---|
| [0d9546] | 301 |   mol = NULL;
 | 
|---|
 | 302 | }
 | 
|---|
 | 303 | 
 | 
|---|
| [e41c48] | 304 | molecule* atom::getMolecule() const {
 | 
|---|
| [c084cc] | 305 |   return mol;
 | 
|---|
 | 306 | }
 | 
|---|
 | 307 | 
 | 
|---|
| [6cfa36] | 308 | void atom::removeFromMolecule(){
 | 
|---|
 | 309 |   if(mol){
 | 
|---|
| [c32d21] | 310 |     if(mol->containsAtom(this)){
 | 
|---|
 | 311 |       signOff(mol, AtomObservable::PositionChanged);
 | 
|---|
| [6cfa36] | 312 |       mol->erase(this);
 | 
|---|
 | 313 |     }
 | 
|---|
 | 314 |     mol=0;
 | 
|---|
 | 315 |   }
 | 
|---|
| [1f8337] | 316 | }
 | 
|---|
 | 317 | 
 | 
|---|
| [560bbe] | 318 | bool atom::changeNr(const int newNr)
 | 
|---|
 | 319 | {
 | 
|---|
 | 320 |   if ((mol) && (mol->changeAtomNr(getNr(),newNr,this))) {
 | 
|---|
 | 321 |     return true;
 | 
|---|
 | 322 |   } else{
 | 
|---|
 | 323 |     return false;
 | 
|---|
 | 324 |   }
 | 
|---|
 | 325 | }
 | 
|---|
 | 326 | 
 | 
|---|
| [e8a21f] | 327 | int atom::getNr() const{
 | 
|---|
| [735b1c] | 328 |   return ParticleInfo::getNr();
 | 
|---|
| [e8a21f] | 329 | }
 | 
|---|
| [6cfa36] | 330 | 
 | 
|---|
| [88d586] | 331 | atom* NewAtom(atomId_t _id){
 | 
|---|
| [11f0fa] | 332 |   atom * res = new atom();
 | 
|---|
| [88d586] | 333 |   res->setId(_id);
 | 
|---|
 | 334 |   return res;
 | 
|---|
| [46d958] | 335 | }
 | 
|---|
 | 336 | 
 | 
|---|
| [88d586] | 337 | void DeleteAtom(atom* atom){
 | 
|---|
| [46d958] | 338 |   delete atom;
 | 
|---|
| [e5f64de] | 339 | }
 | 
|---|
 | 340 | 
 | 
|---|
 | 341 | bool compareAtomElements(atom* atom1,atom* atom2){
 | 
|---|
| [ed26ae] | 342 |   return atom1->getType()->getAtomicNumber() < atom2->getType()->getAtomicNumber();
 | 
|---|
| [46d958] | 343 | }
 | 
|---|
| [d05088] | 344 | /*
 | 
|---|
 | 345 | void atom::update(Observable *publisher)
 | 
|---|
 | 346 | {}
 | 
|---|
 | 347 | 
 | 
|---|
 | 348 | void atom::recieveNotification(Observable *publisher, Notification_ptr notification)
 | 
|---|
 | 349 | {
 | 
|---|
 | 350 |   ASSERT(0, "atom::recieveNotification() - we are not signed on to any notifications.");
 | 
|---|
 | 351 | }
 | 
|---|
 | 352 | */
 | 
|---|
 | 353 | void atom::subjectKilled(Observable *publisher)
 | 
|---|
 | 354 | {
 | 
|---|
 | 355 |   // as publisher has been half-deallocated (Observable is one of the base classes, hence
 | 
|---|
 | 356 |   // becomes destroyed latest), we cannot senibly cast it anymore.
 | 
|---|
 | 357 |   // Hence, we simply have to check here whether it is NOT one of the other instances
 | 
|---|
 | 358 |   // we are signed on to.
 | 
|---|
 | 359 |   father = this;
 | 
|---|
 | 360 |   // no need to sign off
 | 
|---|
 | 361 | }
 | 
|---|