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