| [bcf653] | 1 | /* | 
|---|
|  | 2 | * Project: MoleCuilder | 
|---|
|  | 3 | * Description: creates and alters molecular systems | 
|---|
| [0aa122] | 4 | * Copyright (C)  2010-2012 University of Bonn. All rights reserved. | 
|---|
| [94d5ac6] | 5 | * | 
|---|
|  | 6 | * | 
|---|
|  | 7 | *   This file is part of MoleCuilder. | 
|---|
|  | 8 | * | 
|---|
|  | 9 | *    MoleCuilder is free software: you can redistribute it and/or modify | 
|---|
|  | 10 | *    it under the terms of the GNU General Public License as published by | 
|---|
|  | 11 | *    the Free Software Foundation, either version 2 of the License, or | 
|---|
|  | 12 | *    (at your option) any later version. | 
|---|
|  | 13 | * | 
|---|
|  | 14 | *    MoleCuilder is distributed in the hope that it will be useful, | 
|---|
|  | 15 | *    but WITHOUT ANY WARRANTY; without even the implied warranty of | 
|---|
|  | 16 | *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the | 
|---|
|  | 17 | *    GNU General Public License for more details. | 
|---|
|  | 18 | * | 
|---|
|  | 19 | *    You should have received a copy of the GNU General Public License | 
|---|
|  | 20 | *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>. | 
|---|
| [bcf653] | 21 | */ | 
|---|
|  | 22 |  | 
|---|
| [14de469] | 23 | /** \file atom.cpp | 
|---|
| [1907a7] | 24 | * | 
|---|
| [14de469] | 25 | * Function implementations for the class atom. | 
|---|
| [1907a7] | 26 | * | 
|---|
| [14de469] | 27 | */ | 
|---|
|  | 28 |  | 
|---|
| [bf3817] | 29 | // include config.h | 
|---|
|  | 30 | #ifdef HAVE_CONFIG_H | 
|---|
|  | 31 | #include <config.h> | 
|---|
|  | 32 | #endif | 
|---|
|  | 33 |  | 
|---|
| [ad011c] | 34 | #include "CodePatterns/MemDebug.hpp" | 
|---|
| [112b09] | 35 |  | 
|---|
| [357fba] | 36 | #include "atom.hpp" | 
|---|
| [129204] | 37 | #include "Bond/bond.hpp" | 
|---|
| [e2373df] | 38 | #include "CodePatterns/Log.hpp" | 
|---|
| [4a7776a] | 39 | #include "config.hpp" | 
|---|
| [3bdb6d] | 40 | #include "Element/element.hpp" | 
|---|
| [57f243] | 41 | #include "LinearAlgebra/Vector.hpp" | 
|---|
| [d346b6] | 42 | #include "World.hpp" | 
|---|
| [6cfa36] | 43 | #include "molecule.hpp" | 
|---|
| [c550dd] | 44 | #include "Shapes/Shape.hpp" | 
|---|
| [a0064e] | 45 |  | 
|---|
| [36166d] | 46 | #include <iomanip> | 
|---|
| [0ba410] | 47 | #include <iostream> | 
|---|
| [36166d] | 48 |  | 
|---|
| [14de469] | 49 | /************************************* Functions for class atom *************************************/ | 
|---|
|  | 50 |  | 
|---|
| [70ff32] | 51 |  | 
|---|
| [46d958] | 52 | atom::atom() : | 
|---|
| [97b825] | 53 | father(this), | 
|---|
| [5309ba] | 54 | sort(&Nr), | 
|---|
| [97b825] | 55 | mol(0) | 
|---|
| [d74077] | 56 | {}; | 
|---|
| [14de469] | 57 |  | 
|---|
| [ea0c8b] | 58 | atom::atom(atom *pointer) : | 
|---|
|  | 59 | ParticleInfo(*pointer), | 
|---|
| [c742fe] | 60 | AtomInfo(*pointer), | 
|---|
| [97b825] | 61 | father(pointer), | 
|---|
| [5309ba] | 62 | sort(&Nr), | 
|---|
| [9df680] | 63 | mol(0) | 
|---|
| [2319ed] | 64 | { | 
|---|
| [443547] | 65 | AtomicPosition = pointer->AtomicPosition; // copy trajectory of coordination | 
|---|
|  | 66 | AtomicVelocity = pointer->AtomicVelocity; // copy trajectory of velocity | 
|---|
|  | 67 | AtomicForce = pointer->AtomicForce; | 
|---|
| [b453f9] | 68 | }; | 
|---|
| [2319ed] | 69 |  | 
|---|
| [46d958] | 70 | atom *atom::clone(){ | 
|---|
| [68f03d] | 71 | atom *res = new atom(this); | 
|---|
| [23b547] | 72 | World::getInstance().registerAtom(res); | 
|---|
| [46d958] | 73 | return res; | 
|---|
|  | 74 | } | 
|---|
|  | 75 |  | 
|---|
| [2319ed] | 76 |  | 
|---|
| [14de469] | 77 | /** Destructor of class atom. | 
|---|
|  | 78 | */ | 
|---|
| [1907a7] | 79 | atom::~atom() | 
|---|
| [14de469] | 80 | { | 
|---|
| [6cfa36] | 81 | removeFromMolecule(); | 
|---|
| [14de469] | 82 | }; | 
|---|
| [e2373df] | 83 |  | 
|---|
|  | 84 |  | 
|---|
|  | 85 | void atom::UpdateSteps() | 
|---|
|  | 86 | { | 
|---|
|  | 87 | LOG(4,"atom::UpdateSteps() called."); | 
|---|
|  | 88 | // append to position, velocity and force vector | 
|---|
|  | 89 | AtomInfo::AppendTrajectoryStep(); | 
|---|
| [1e6249] | 90 | // append to ListOfBonds vector | 
|---|
|  | 91 | BondedParticleInfo::AppendTrajectoryStep(); | 
|---|
| [e2373df] | 92 | } | 
|---|
|  | 93 |  | 
|---|
| [59fff1] | 94 | atom *atom::GetTrueFather() | 
|---|
|  | 95 | { | 
|---|
|  | 96 | const atom *father = const_cast<const atom *>(this)->GetTrueFather(); | 
|---|
|  | 97 | return const_cast<atom *>(father); | 
|---|
|  | 98 | } | 
|---|
|  | 99 |  | 
|---|
|  | 100 | const atom *atom::GetTrueFather() const | 
|---|
| [14de469] | 101 | { | 
|---|
| [215df0] | 102 | if(father == this){ // top most father is the one that points on itself | 
|---|
|  | 103 | return this; | 
|---|
|  | 104 | } | 
|---|
|  | 105 | else if(!father) { | 
|---|
|  | 106 | return 0; | 
|---|
|  | 107 | } | 
|---|
|  | 108 | else { | 
|---|
|  | 109 | return father->GetTrueFather(); | 
|---|
|  | 110 | } | 
|---|
| [14de469] | 111 | }; | 
|---|
|  | 112 |  | 
|---|
| [e65246] | 113 | /** Sets father to itself or its father in case of copying a molecule. | 
|---|
|  | 114 | */ | 
|---|
|  | 115 | void atom::CorrectFather() | 
|---|
|  | 116 | { | 
|---|
| [2e352f] | 117 | if (father->father != father)   // same atom in copy's father points to itself | 
|---|
|  | 118 | //    father = this;  // set father to itself (copy of a whole molecule) | 
|---|
|  | 119 | //  else | 
|---|
| [e65246] | 120 | father = father->father;  // set father to original's father | 
|---|
|  | 121 |  | 
|---|
|  | 122 | }; | 
|---|
|  | 123 |  | 
|---|
| [b453f9] | 124 | void atom::EqualsFather ( const atom *ptr, const atom **res ) const | 
|---|
| [e65246] | 125 | { | 
|---|
|  | 126 | if ( ptr == father ) | 
|---|
|  | 127 | *res = this; | 
|---|
|  | 128 | }; | 
|---|
|  | 129 |  | 
|---|
| [00abfc] | 130 | bool atom::isFather(const atom *ptr){ | 
|---|
|  | 131 | return ptr==father; | 
|---|
|  | 132 | } | 
|---|
|  | 133 |  | 
|---|
| [e138de] | 134 | bool atom::OutputIndexed(ofstream * const out, const int ElementNo, const int AtomNo, const char *comment) const | 
|---|
| [14de469] | 135 | { | 
|---|
|  | 136 | if (out != NULL) { | 
|---|
|  | 137 | *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t"  << fixed << setprecision(9) << showpoint; | 
|---|
| [d74077] | 138 | *out << at(0) << "\t" << at(1) << "\t" << at(2); | 
|---|
| [6625c3] | 139 | *out << "\t" << (int)(getFixedIon()); | 
|---|
| [bce72c] | 140 | if (getAtomicVelocity().Norm() > MYEPSILON) | 
|---|
|  | 141 | *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t"; | 
|---|
| [437922] | 142 | if (comment != NULL) | 
|---|
|  | 143 | *out << " # " << comment << endl; | 
|---|
| [e9f8f9] | 144 | else | 
|---|
| [735b1c] | 145 | *out << " # molecule nr " << getNr() << endl; | 
|---|
| [e9f8f9] | 146 | return true; | 
|---|
|  | 147 | } else | 
|---|
|  | 148 | return false; | 
|---|
|  | 149 | }; | 
|---|
| [b453f9] | 150 |  | 
|---|
| [0ba410] | 151 | bool atom::OutputArrayIndexed(ostream * const out,const enumeration<const element*> &elementLookup, int *AtomNo, const char *comment) const | 
|---|
| [e9f8f9] | 152 | { | 
|---|
| [83f176] | 153 | AtomNo[getType()->getAtomicNumber()]++;  // increment number | 
|---|
| [e9f8f9] | 154 | if (out != NULL) { | 
|---|
| [8f4df1] | 155 | const element *elemental = getType(); | 
|---|
|  | 156 | ASSERT(elementLookup.there.find(elemental)!=elementLookup.there.end(),"Type of this atom was not in the formula upon enumeration"); | 
|---|
| [83f176] | 157 | *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[elemental->getAtomicNumber()] << "\t"  << fixed << setprecision(9) << showpoint; | 
|---|
| [d74077] | 158 | *out << at(0) << "\t" << at(1) << "\t" << at(2); | 
|---|
| [6625c3] | 159 | *out << "\t" << getFixedIon(); | 
|---|
| [bce72c] | 160 | if (getAtomicVelocity().Norm() > MYEPSILON) | 
|---|
|  | 161 | *out << "\t" << scientific << setprecision(6) << getAtomicVelocity()[0] << "\t" << getAtomicVelocity()[1] << "\t" << getAtomicVelocity()[2] << "\t"; | 
|---|
| [e9f8f9] | 162 | if (comment != NULL) | 
|---|
|  | 163 | *out << " # " << comment << endl; | 
|---|
| [437922] | 164 | else | 
|---|
| [735b1c] | 165 | *out << " # molecule nr " << getNr() << endl; | 
|---|
| [14de469] | 166 | return true; | 
|---|
|  | 167 | } else | 
|---|
|  | 168 | return false; | 
|---|
|  | 169 | }; | 
|---|
|  | 170 |  | 
|---|
| [b453f9] | 171 | bool atom::Compare(const atom &ptr) const | 
|---|
| [4455f4] | 172 | { | 
|---|
| [735b1c] | 173 | if (getNr() < ptr.getNr()) | 
|---|
| [4455f4] | 174 | return true; | 
|---|
|  | 175 | else | 
|---|
|  | 176 | return false; | 
|---|
|  | 177 | }; | 
|---|
|  | 178 |  | 
|---|
| [b453f9] | 179 | double atom::DistanceSquaredToVector(const Vector &origin) const | 
|---|
| [4455f4] | 180 | { | 
|---|
| [d74077] | 181 | return DistanceSquared(origin); | 
|---|
| [4455f4] | 182 | }; | 
|---|
|  | 183 |  | 
|---|
| [b453f9] | 184 | double atom::DistanceToVector(const Vector &origin) const | 
|---|
| [4455f4] | 185 | { | 
|---|
| [d74077] | 186 | return distance(origin); | 
|---|
| [4455f4] | 187 | }; | 
|---|
|  | 188 |  | 
|---|
|  | 189 | void atom::InitComponentNr() | 
|---|
|  | 190 | { | 
|---|
|  | 191 | if (ComponentNr != NULL) | 
|---|
| [920c70] | 192 | delete[](ComponentNr); | 
|---|
| [9d83b6] | 193 | const BondList& ListOfBonds = getListOfBonds(); | 
|---|
| [920c70] | 194 | ComponentNr = new int[ListOfBonds.size()+1]; | 
|---|
| [4455f4] | 195 | for (int i=ListOfBonds.size()+1;i--;) | 
|---|
|  | 196 | ComponentNr[i] = -1; | 
|---|
| [14b65e] | 197 | }; | 
|---|
|  | 198 |  | 
|---|
|  | 199 | void atom::resetGraphNr(){ | 
|---|
|  | 200 | GraphNr=-1; | 
|---|
|  | 201 | } | 
|---|
| [4455f4] | 202 |  | 
|---|
| [d74077] | 203 | std::ostream & atom::operator << (std::ostream &ost) const | 
|---|
|  | 204 | { | 
|---|
|  | 205 | ParticleInfo::operator<<(ost); | 
|---|
|  | 206 | ost << "," << getPosition(); | 
|---|
|  | 207 | return ost; | 
|---|
|  | 208 | } | 
|---|
|  | 209 |  | 
|---|
|  | 210 | std::ostream & operator << (std::ostream &ost, const atom &a) | 
|---|
|  | 211 | { | 
|---|
|  | 212 | a.ParticleInfo::operator<<(ost); | 
|---|
|  | 213 | ost << "," << a.getPosition(); | 
|---|
|  | 214 | return ost; | 
|---|
|  | 215 | } | 
|---|
| [4455f4] | 216 |  | 
|---|
|  | 217 | bool operator < (atom &a, atom &b) | 
|---|
|  | 218 | { | 
|---|
|  | 219 | return a.Compare(b); | 
|---|
|  | 220 | }; | 
|---|
|  | 221 |  | 
|---|
| [46d958] | 222 | World *atom::getWorld(){ | 
|---|
|  | 223 | return world; | 
|---|
|  | 224 | } | 
|---|
|  | 225 |  | 
|---|
|  | 226 | void atom::setWorld(World* _world){ | 
|---|
|  | 227 | world = _world; | 
|---|
|  | 228 | } | 
|---|
|  | 229 |  | 
|---|
| [88d586] | 230 | bool atom::changeId(atomId_t newId){ | 
|---|
|  | 231 | // first we move ourselves in the world | 
|---|
|  | 232 | // the world lets us know if that succeeded | 
|---|
| [4f7f0bf] | 233 | if(world->changeAtomId(id,newId,this)){ | 
|---|
|  | 234 | OBSERVE; | 
|---|
|  | 235 | id = newId; | 
|---|
|  | 236 | NOTIFY(IndexChanged); | 
|---|
| [88d586] | 237 | return true; | 
|---|
|  | 238 | } | 
|---|
|  | 239 | else{ | 
|---|
|  | 240 | return false; | 
|---|
|  | 241 | } | 
|---|
|  | 242 | } | 
|---|
|  | 243 |  | 
|---|
|  | 244 | void atom::setId(atomId_t _id) { | 
|---|
| [46d958] | 245 | id=_id; | 
|---|
|  | 246 | } | 
|---|
|  | 247 |  | 
|---|
| [ad2b411] | 248 | atomId_t atom::getId() const { | 
|---|
| [46d958] | 249 | return id; | 
|---|
|  | 250 | } | 
|---|
|  | 251 |  | 
|---|
| [6cfa36] | 252 | void atom::setMolecule(molecule *_mol){ | 
|---|
|  | 253 | // take this atom from the old molecule | 
|---|
|  | 254 | removeFromMolecule(); | 
|---|
| [3867a7] | 255 | mol = _mol; | 
|---|
|  | 256 | if ((mol) && (!mol->containsAtom(this))) | 
|---|
| [dddbfe] | 257 | mol->insert(this); | 
|---|
| [6cfa36] | 258 | } | 
|---|
|  | 259 |  | 
|---|
| [0d9546] | 260 | void atom::unsetMolecule() | 
|---|
|  | 261 | { | 
|---|
|  | 262 | // take this atom from the old molecule | 
|---|
|  | 263 | ASSERT(!mol->containsAtom(this), | 
|---|
|  | 264 | "atom::unsetMolecule() - old molecule "+toString(mol)+" still contains us!"); | 
|---|
|  | 265 | mol = NULL; | 
|---|
|  | 266 | } | 
|---|
|  | 267 |  | 
|---|
| [e41c48] | 268 | molecule* atom::getMolecule() const { | 
|---|
| [c084cc] | 269 | return mol; | 
|---|
|  | 270 | } | 
|---|
|  | 271 |  | 
|---|
| [6cfa36] | 272 | void atom::removeFromMolecule(){ | 
|---|
|  | 273 | if(mol){ | 
|---|
|  | 274 | if(mol->containsAtom(this)){ | 
|---|
|  | 275 | mol->erase(this); | 
|---|
|  | 276 | } | 
|---|
|  | 277 | mol=0; | 
|---|
|  | 278 | } | 
|---|
| [1f8337] | 279 | } | 
|---|
|  | 280 |  | 
|---|
| [560bbe] | 281 | bool atom::changeNr(const int newNr) | 
|---|
|  | 282 | { | 
|---|
|  | 283 | if ((mol) && (mol->changeAtomNr(getNr(),newNr,this))) { | 
|---|
|  | 284 | return true; | 
|---|
|  | 285 | } else{ | 
|---|
|  | 286 | return false; | 
|---|
|  | 287 | } | 
|---|
|  | 288 | } | 
|---|
|  | 289 |  | 
|---|
| [e8a21f] | 290 | int atom::getNr() const{ | 
|---|
| [735b1c] | 291 | return ParticleInfo::getNr(); | 
|---|
| [e8a21f] | 292 | } | 
|---|
| [6cfa36] | 293 |  | 
|---|
| [88d586] | 294 | atom* NewAtom(atomId_t _id){ | 
|---|
|  | 295 | atom * res =new atom(); | 
|---|
|  | 296 | res->setId(_id); | 
|---|
|  | 297 | return res; | 
|---|
| [46d958] | 298 | } | 
|---|
|  | 299 |  | 
|---|
| [88d586] | 300 | void DeleteAtom(atom* atom){ | 
|---|
| [46d958] | 301 | delete atom; | 
|---|
| [e5f64de] | 302 | } | 
|---|
|  | 303 |  | 
|---|
|  | 304 | bool compareAtomElements(atom* atom1,atom* atom2){ | 
|---|
| [ed26ae] | 305 | return atom1->getType()->getAtomicNumber() < atom2->getType()->getAtomicNumber(); | 
|---|
| [46d958] | 306 | } | 
|---|