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) 2014 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 | /*
|
---|
25 | * atom_atominfo.cpp
|
---|
26 | *
|
---|
27 | * Created on: Oct 19, 2009
|
---|
28 | * Author: heber
|
---|
29 | */
|
---|
30 |
|
---|
31 | // include config.h
|
---|
32 | #ifdef HAVE_CONFIG_H
|
---|
33 | #include <config.h>
|
---|
34 | #endif
|
---|
35 |
|
---|
36 | //#include "CodePatterns/MemDebug.hpp"
|
---|
37 |
|
---|
38 | #include "CodePatterns/Verbose.hpp"
|
---|
39 |
|
---|
40 | #include "atom_atominfo.hpp"
|
---|
41 | #include "CodePatterns/Log.hpp"
|
---|
42 | #include "config.hpp"
|
---|
43 | #include "Element/element.hpp"
|
---|
44 | #include "Element/periodentafel.hpp"
|
---|
45 | #include "Fragmentation/ForceMatrix.hpp"
|
---|
46 | #include "World.hpp"
|
---|
47 | #include "WorldTime.hpp"
|
---|
48 |
|
---|
49 | #include <iomanip>
|
---|
50 |
|
---|
51 | /** Constructor of class AtomInfo.
|
---|
52 | */
|
---|
53 | AtomInfo::AtomInfo() :
|
---|
54 | AtomicElement(1),
|
---|
55 | FixedIon(false),
|
---|
56 | charge(0.)
|
---|
57 | {
|
---|
58 | AtomicPosition.insert( std::make_pair(0, zeroVec) );
|
---|
59 | AtomicVelocity.insert( std::make_pair(0, zeroVec) );
|
---|
60 | AtomicForce.insert( std::make_pair(0, zeroVec) );
|
---|
61 | }
|
---|
62 |
|
---|
63 | /** Copy constructor of class AtomInfo.
|
---|
64 | */
|
---|
65 | AtomInfo::AtomInfo(const AtomInfo &_atom) :
|
---|
66 | AtomicPosition(_atom.AtomicPosition),
|
---|
67 | AtomicVelocity(_atom.AtomicVelocity),
|
---|
68 | AtomicForce(_atom.AtomicForce),
|
---|
69 | AtomicElement(_atom.AtomicElement),
|
---|
70 | FixedIon(_atom.FixedIon),
|
---|
71 | charge(_atom.charge)
|
---|
72 | {
|
---|
73 | }
|
---|
74 |
|
---|
75 | AtomInfo::AtomInfo(const VectorInterface &_v) :
|
---|
76 | AtomicElement(1),
|
---|
77 | FixedIon(false),
|
---|
78 | charge(0.)
|
---|
79 | {
|
---|
80 | AtomicPosition.insert( std::make_pair(0, _v.getPosition()) );
|
---|
81 | AtomicVelocity.insert( std::make_pair(0, zeroVec) );
|
---|
82 | AtomicForce.insert( std::make_pair(0, zeroVec) );
|
---|
83 | };
|
---|
84 |
|
---|
85 | /** Destructor of class AtomInfo.
|
---|
86 | */
|
---|
87 | AtomInfo::~AtomInfo()
|
---|
88 | {
|
---|
89 | };
|
---|
90 |
|
---|
91 | void AtomInfo::AppendTrajectoryStep(const unsigned int _step)
|
---|
92 | {
|
---|
93 | NOTIFY(TrajectoryChanged);
|
---|
94 | AtomicPosition.insert( std::make_pair(_step, zeroVec) );
|
---|
95 | AtomicVelocity.insert( std::make_pair(_step, zeroVec) );
|
---|
96 | AtomicForce.insert( std::make_pair(_step, zeroVec) );
|
---|
97 | LOG(5,"AtomInfo::AppendTrajectoryStep() called, size is ("
|
---|
98 | << AtomicPosition.size() << ","
|
---|
99 | << AtomicVelocity.size() << ","
|
---|
100 | << AtomicForce.size() << ")");
|
---|
101 | }
|
---|
102 |
|
---|
103 | void AtomInfo::removeTrajectoryStep(const unsigned int _step)
|
---|
104 | {
|
---|
105 | NOTIFY(TrajectoryChanged);
|
---|
106 | AtomicPosition.erase(_step);
|
---|
107 | AtomicVelocity.erase(_step);
|
---|
108 | AtomicForce.erase(_step);
|
---|
109 | LOG(5,"AtomInfo::removeTrajectoryStep() called, size is ("
|
---|
110 | << AtomicPosition.size() << ","
|
---|
111 | << AtomicVelocity.size() << ","
|
---|
112 | << AtomicForce.size() << ")");
|
---|
113 | }
|
---|
114 |
|
---|
115 | const element *AtomInfo::getType() const
|
---|
116 | {
|
---|
117 | const element *elem = World::getInstance().getPeriode()->FindElement(AtomicElement);
|
---|
118 | return elem;
|
---|
119 | }
|
---|
120 |
|
---|
121 | const element &AtomInfo::getElement() const
|
---|
122 | {
|
---|
123 | const element &elem = *World::getInstance().getPeriode()->FindElement(AtomicElement);
|
---|
124 | return elem;
|
---|
125 | }
|
---|
126 |
|
---|
127 | atomicNumber_t AtomInfo::getElementNo() const
|
---|
128 | {
|
---|
129 | return AtomicElement;
|
---|
130 | }
|
---|
131 |
|
---|
132 | const std::string &AtomInfo::getParticleName() const
|
---|
133 | {
|
---|
134 | return particlename;
|
---|
135 | }
|
---|
136 |
|
---|
137 | void AtomInfo::setParticleName(const std::string & _name)
|
---|
138 | {
|
---|
139 | particlename = _name;
|
---|
140 | }
|
---|
141 |
|
---|
142 | const double& AtomInfo::operator[](size_t i) const
|
---|
143 | {
|
---|
144 | return atStep(i, WorldTime::getTime());
|
---|
145 | }
|
---|
146 |
|
---|
147 | const double& AtomInfo::at(size_t i) const
|
---|
148 | {
|
---|
149 | return atStep(i, WorldTime::getTime());
|
---|
150 | }
|
---|
151 |
|
---|
152 | const double& AtomInfo::atStep(size_t i, unsigned int _step) const
|
---|
153 | {
|
---|
154 | ASSERT(!AtomicPosition.empty(),
|
---|
155 | "AtomInfo::operator[]() - AtomicPosition is empty.");
|
---|
156 | VectorTrajectory_t::const_iterator iter =
|
---|
157 | AtomicPosition.lower_bound(_step);
|
---|
158 | return iter->second[i];
|
---|
159 | }
|
---|
160 |
|
---|
161 | void AtomInfo::set(size_t i, const double value)
|
---|
162 | {
|
---|
163 | OBSERVE;
|
---|
164 | NOTIFY(AtomObservable::PositionChanged);
|
---|
165 | VectorTrajectory_t::iterator iter = AtomicPosition.find(WorldTime::getTime());
|
---|
166 | if (iter != AtomicPosition.end()) {
|
---|
167 | iter->second[i] = value;
|
---|
168 | } else {
|
---|
169 | Vector newPos;
|
---|
170 | newPos[i] = value;
|
---|
171 | #ifndef NDEBUG
|
---|
172 | std::pair<VectorTrajectory_t::iterator, bool> inserter =
|
---|
173 | #endif
|
---|
174 | AtomicPosition.insert( std::make_pair(WorldTime::getTime(), newPos) );
|
---|
175 | ASSERT( inserter.second,
|
---|
176 | "AtomInfo::set() - time step "+toString(WorldTime::getTime())
|
---|
177 | +" present after all?");
|
---|
178 | }
|
---|
179 | }
|
---|
180 |
|
---|
181 | void AtomInfo::setAtStep(size_t i, unsigned int _step, const double value)
|
---|
182 | {
|
---|
183 | OBSERVE;
|
---|
184 | NOTIFY(AtomObservable::PositionChanged);
|
---|
185 | VectorTrajectory_t::iterator iter = AtomicPosition.find(_step);
|
---|
186 | if (iter != AtomicPosition.end()) {
|
---|
187 | iter->second[i] = value;
|
---|
188 | } else {
|
---|
189 | Vector newPos;
|
---|
190 | newPos[i] = value;
|
---|
191 | #ifndef NDEBUG
|
---|
192 | std::pair<VectorTrajectory_t::iterator, bool> inserter =
|
---|
193 | #endif
|
---|
194 | AtomicPosition.insert( std::make_pair(_step, newPos) );
|
---|
195 | ASSERT( inserter.second,
|
---|
196 | "AtomInfo::setAtStep() - time step "+toString(_step)
|
---|
197 | +" present after all?");
|
---|
198 | }
|
---|
199 | }
|
---|
200 |
|
---|
201 | /** Helps to determine whether the current step really exists or getPosition() has just
|
---|
202 | * delivered the one closest to it in the past.
|
---|
203 | *
|
---|
204 | * \param _step step to check
|
---|
205 | * \param true - step exists, false - step does not exist, getPosition() delivers closest
|
---|
206 | */
|
---|
207 | bool AtomInfo::isStepPresent(const unsigned int _step) const
|
---|
208 | {
|
---|
209 | VectorTrajectory_t::const_iterator iter =
|
---|
210 | AtomicPosition.find(_step);
|
---|
211 | return iter != AtomicPosition.end();
|
---|
212 | }
|
---|
213 |
|
---|
214 | const Vector& AtomInfo::getPosition() const
|
---|
215 | {
|
---|
216 | return getPositionAtStep(WorldTime::getTime());
|
---|
217 | }
|
---|
218 |
|
---|
219 | const Vector& AtomInfo::getPositionAtStep(const unsigned int _step) const
|
---|
220 | {
|
---|
221 | ASSERT(!AtomicPosition.empty(),
|
---|
222 | "AtomInfo::operator[]() - AtomicPosition is empty.");
|
---|
223 | VectorTrajectory_t::const_iterator iter =
|
---|
224 | AtomicPosition.lower_bound(_step);
|
---|
225 | return iter->second;
|
---|
226 | }
|
---|
227 |
|
---|
228 | void AtomInfo::setType(const element* _type)
|
---|
229 | {
|
---|
230 | OBSERVE;
|
---|
231 | NOTIFY(AtomObservable::ElementChanged);
|
---|
232 | AtomicElement = _type->getAtomicNumber();
|
---|
233 | }
|
---|
234 |
|
---|
235 | void AtomInfo::setType(const int Z)
|
---|
236 | {
|
---|
237 | const element *elem = World::getInstance().getPeriode()->FindElement(Z);
|
---|
238 | setType(elem);
|
---|
239 | }
|
---|
240 |
|
---|
241 | const Vector& AtomInfo::getAtomicVelocity() const
|
---|
242 | {
|
---|
243 | return getAtomicVelocityAtStep(WorldTime::getTime());
|
---|
244 | }
|
---|
245 |
|
---|
246 | const Vector& AtomInfo::getAtomicVelocityAtStep(const unsigned int _step) const
|
---|
247 | {
|
---|
248 | ASSERT(!AtomicVelocity.empty(),
|
---|
249 | "AtomInfo::operator[]() - AtomicVelocity is empty.");
|
---|
250 | VectorTrajectory_t::const_iterator iter =
|
---|
251 | AtomicVelocity.lower_bound(_step);
|
---|
252 | // special, we only interpolate between present time steps not into the future
|
---|
253 | if (_step > AtomicVelocity.begin()->first)
|
---|
254 | return zeroVec;
|
---|
255 | else
|
---|
256 | return iter->second;
|
---|
257 | }
|
---|
258 |
|
---|
259 | void AtomInfo::setAtomicVelocity(const Vector &_newvelocity)
|
---|
260 | {
|
---|
261 | setAtomicVelocityAtStep(WorldTime::getTime(), _newvelocity);
|
---|
262 | }
|
---|
263 |
|
---|
264 | void AtomInfo::setAtomicVelocityAtStep(const unsigned int _step, const Vector &_newvelocity)
|
---|
265 | {
|
---|
266 | OBSERVE;
|
---|
267 | VectorTrajectory_t::iterator iter = AtomicVelocity.find(_step);
|
---|
268 | if (iter != AtomicVelocity.end()) {
|
---|
269 | iter->second = _newvelocity;
|
---|
270 | } else {
|
---|
271 | #ifndef NDEBUG
|
---|
272 | std::pair<VectorTrajectory_t::iterator, bool> inserter =
|
---|
273 | #endif
|
---|
274 | AtomicVelocity.insert( std::make_pair(_step, _newvelocity) );
|
---|
275 | ASSERT( inserter.second,
|
---|
276 | "AtomInfo::set() - time step "+toString(_step)
|
---|
277 | +" present after all?");
|
---|
278 | }
|
---|
279 | if (WorldTime::getTime() == _step)
|
---|
280 | NOTIFY(AtomObservable::VelocityChanged);
|
---|
281 | }
|
---|
282 |
|
---|
283 | const Vector& AtomInfo::getAtomicForce() const
|
---|
284 | {
|
---|
285 | return getAtomicForceAtStep(WorldTime::getTime());
|
---|
286 | }
|
---|
287 |
|
---|
288 | const Vector& AtomInfo::getAtomicForceAtStep(const unsigned int _step) const
|
---|
289 | {
|
---|
290 | ASSERT(!AtomicForce.empty(),
|
---|
291 | "AtomInfo::operator[]() - AtomicForce is empty.");
|
---|
292 | VectorTrajectory_t::const_iterator iter =
|
---|
293 | AtomicForce.lower_bound(_step);
|
---|
294 | // special, we only interpolate between present time steps not into the future
|
---|
295 | if (_step > AtomicForce.begin()->first)
|
---|
296 | return zeroVec;
|
---|
297 | else
|
---|
298 | return iter->second;
|
---|
299 | }
|
---|
300 |
|
---|
301 | void AtomInfo::setAtomicForce(const Vector &_newforce)
|
---|
302 | {
|
---|
303 | setAtomicForceAtStep(WorldTime::getTime(), _newforce);
|
---|
304 | }
|
---|
305 |
|
---|
306 | void AtomInfo::setAtomicForceAtStep(const unsigned int _step, const Vector &_newforce)
|
---|
307 | {
|
---|
308 | OBSERVE;
|
---|
309 | VectorTrajectory_t::iterator iter = AtomicForce.find(_step);
|
---|
310 | if (iter != AtomicForce.end()) {
|
---|
311 | iter->second = _newforce;
|
---|
312 | } else {
|
---|
313 | #ifndef NDEBUG
|
---|
314 | std::pair<VectorTrajectory_t::iterator, bool> inserter =
|
---|
315 | #endif
|
---|
316 | AtomicForce.insert( std::make_pair(_step, _newforce) );
|
---|
317 | ASSERT( inserter.second,
|
---|
318 | "AtomInfo::set() - time step "+toString(_step)
|
---|
319 | +" present after all?");
|
---|
320 | }
|
---|
321 | if (WorldTime::getTime() == _step)
|
---|
322 | NOTIFY(AtomObservable::ForceChanged);
|
---|
323 | }
|
---|
324 |
|
---|
325 | bool AtomInfo::getFixedIon() const
|
---|
326 | {
|
---|
327 | return FixedIon;
|
---|
328 | }
|
---|
329 |
|
---|
330 | void AtomInfo::setFixedIon(const bool _fixedion)
|
---|
331 | {
|
---|
332 | OBSERVE;
|
---|
333 | NOTIFY(AtomObservable::PropertyChanged);
|
---|
334 | FixedIon = _fixedion;
|
---|
335 | }
|
---|
336 |
|
---|
337 | void AtomInfo::setPosition(const Vector& _vector)
|
---|
338 | {
|
---|
339 | setPositionAtStep(WorldTime::getTime(), _vector);
|
---|
340 | }
|
---|
341 |
|
---|
342 | void AtomInfo::setPositionAtStep(unsigned int _step, const Vector& _vector)
|
---|
343 | {
|
---|
344 | OBSERVE;
|
---|
345 | VectorTrajectory_t::iterator iter = AtomicPosition.find(_step);
|
---|
346 | if (iter != AtomicPosition.end()) {
|
---|
347 | iter->second = _vector;
|
---|
348 | } else {
|
---|
349 | #ifndef NDEBUG
|
---|
350 | std::pair<VectorTrajectory_t::iterator, bool> inserter =
|
---|
351 | #endif
|
---|
352 | AtomicPosition.insert( std::make_pair(_step, _vector) );
|
---|
353 | ASSERT( inserter.second,
|
---|
354 | "AtomInfo::set() - time step "+toString(_step)
|
---|
355 | +" present after all?");
|
---|
356 | }
|
---|
357 | if (WorldTime::getTime() == _step)
|
---|
358 | NOTIFY(AtomObservable::PositionChanged);
|
---|
359 | }
|
---|
360 |
|
---|
361 | const VectorInterface& AtomInfo::operator+=(const Vector& b)
|
---|
362 | {
|
---|
363 | setPosition(getPosition()+b);
|
---|
364 | return *this;
|
---|
365 | }
|
---|
366 |
|
---|
367 | const VectorInterface& AtomInfo::operator-=(const Vector& b)
|
---|
368 | {
|
---|
369 | setPosition(getPosition()-b);
|
---|
370 | return *this;
|
---|
371 | }
|
---|
372 |
|
---|
373 | Vector const AtomInfo::operator+(const Vector& b) const
|
---|
374 | {
|
---|
375 | Vector a(getPosition());
|
---|
376 | a += b;
|
---|
377 | return a;
|
---|
378 | }
|
---|
379 |
|
---|
380 | Vector const AtomInfo::operator-(const Vector& b) const
|
---|
381 | {
|
---|
382 | Vector a(getPosition());
|
---|
383 | a -= b;
|
---|
384 | return a;
|
---|
385 | }
|
---|
386 |
|
---|
387 | double AtomInfo::distance(const Vector &point) const
|
---|
388 | {
|
---|
389 | return getPosition().distance(point);
|
---|
390 | }
|
---|
391 |
|
---|
392 | double AtomInfo::DistanceSquared(const Vector &y) const
|
---|
393 | {
|
---|
394 | return getPosition().DistanceSquared(y);
|
---|
395 | }
|
---|
396 |
|
---|
397 | double AtomInfo::distance(const VectorInterface &_atom) const
|
---|
398 | {
|
---|
399 | return _atom.distance(getPosition());
|
---|
400 | }
|
---|
401 |
|
---|
402 | double AtomInfo::DistanceSquared(const VectorInterface &_atom) const
|
---|
403 | {
|
---|
404 | return _atom.DistanceSquared(getPosition());
|
---|
405 | }
|
---|
406 |
|
---|
407 | VectorInterface &AtomInfo::operator=(const Vector& _vector)
|
---|
408 | {
|
---|
409 | setPosition(_vector);
|
---|
410 | return *this;
|
---|
411 | }
|
---|
412 |
|
---|
413 | void AtomInfo::ScaleAll(const double *factor)
|
---|
414 | {
|
---|
415 | Vector temp(getPosition());
|
---|
416 | temp.ScaleAll(factor);
|
---|
417 | setPosition(temp);
|
---|
418 | }
|
---|
419 |
|
---|
420 | void AtomInfo::ScaleAll(const Vector &factor)
|
---|
421 | {
|
---|
422 | Vector temp(getPosition());
|
---|
423 | temp.ScaleAll(factor);
|
---|
424 | setPosition(temp);
|
---|
425 | }
|
---|
426 |
|
---|
427 | void AtomInfo::Scale(const double factor)
|
---|
428 | {
|
---|
429 | Vector temp(getPosition());
|
---|
430 | temp.Scale(factor);
|
---|
431 | setPosition(temp);
|
---|
432 | }
|
---|
433 |
|
---|
434 | void AtomInfo::Zero()
|
---|
435 | {
|
---|
436 | setPosition(zeroVec);
|
---|
437 | }
|
---|
438 |
|
---|
439 | void AtomInfo::One(const double one)
|
---|
440 | {
|
---|
441 | setPosition(Vector(one,one,one));
|
---|
442 | }
|
---|
443 |
|
---|
444 | void AtomInfo::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
|
---|
445 | {
|
---|
446 | Vector newPos;
|
---|
447 | newPos.LinearCombinationOfVectors(x1,x2,x3,factors);
|
---|
448 | setPosition(newPos);
|
---|
449 | }
|
---|
450 |
|
---|
451 | /**
|
---|
452 | * returns the kinetic energy of this atom at a given time step
|
---|
453 | */
|
---|
454 | double AtomInfo::getKineticEnergy(const unsigned int _step) const
|
---|
455 | {
|
---|
456 | return getMass() * getAtomicVelocityAtStep(_step).NormSquared();
|
---|
457 | }
|
---|
458 |
|
---|
459 | Vector AtomInfo::getMomentum(const unsigned int _step) const
|
---|
460 | {
|
---|
461 | return getMass() * getAtomicVelocityAtStep(_step);
|
---|
462 | }
|
---|
463 |
|
---|
464 | /** Decrease the trajectory if given \a MaxSteps is smaller.
|
---|
465 | * Does nothing if \a MaxSteps is larger than current size.
|
---|
466 | *
|
---|
467 | * \param MaxSteps
|
---|
468 | */
|
---|
469 | void AtomInfo::ResizeTrajectory(size_t MaxSteps)
|
---|
470 | {
|
---|
471 | // mind the reverse ordering due to std::greater, latest time steps are at beginning
|
---|
472 | VectorTrajectory_t::iterator positer = AtomicPosition.lower_bound(MaxSteps);
|
---|
473 | if (positer != AtomicPosition.begin()) {
|
---|
474 | if (positer->first == MaxSteps)
|
---|
475 | --positer;
|
---|
476 | AtomicPosition.erase(AtomicPosition.begin(), positer);
|
---|
477 | }
|
---|
478 | VectorTrajectory_t::iterator veliter = AtomicVelocity.lower_bound(MaxSteps);
|
---|
479 | if (veliter != AtomicVelocity.begin()) {
|
---|
480 | if (veliter->first == MaxSteps)
|
---|
481 | --veliter;
|
---|
482 | AtomicVelocity.erase(AtomicVelocity.begin(), veliter);
|
---|
483 | }
|
---|
484 | VectorTrajectory_t::iterator forceiter = AtomicForce.lower_bound(MaxSteps);
|
---|
485 | if (forceiter != AtomicForce.begin()) {
|
---|
486 | if (forceiter->first == MaxSteps)
|
---|
487 | --forceiter;
|
---|
488 | AtomicForce.erase(AtomicForce.begin(), forceiter);
|
---|
489 | }
|
---|
490 | }
|
---|
491 |
|
---|
492 | size_t AtomInfo::getTrajectorySize() const
|
---|
493 | {
|
---|
494 | // mind greater comp for map here: first element is latest in time steps!
|
---|
495 | return AtomicPosition.begin()->first+1;
|
---|
496 | }
|
---|
497 |
|
---|
498 | double AtomInfo::getMass() const
|
---|
499 | {
|
---|
500 | return getType()->getMass();
|
---|
501 | }
|
---|
502 |
|
---|
503 | /** Helper function to either insert or assign, depending on the element being
|
---|
504 | * present already.
|
---|
505 | *
|
---|
506 | * \param _trajectory vector of Vectors to assign
|
---|
507 | * \param dest step to insert/assign to
|
---|
508 | * \param _newvalue new Vector value
|
---|
509 | */
|
---|
510 | void assignTrajectoryElement(
|
---|
511 | std::map<unsigned int, Vector, std::greater<unsigned int> > &_trajectory,
|
---|
512 | const unsigned int dest,
|
---|
513 | const Vector &_newvalue)
|
---|
514 | {
|
---|
515 | std::pair<std::map<unsigned int, Vector, std::greater<unsigned int> >::iterator, bool> inserter =
|
---|
516 | _trajectory.insert( std::make_pair(dest, _newvalue) );
|
---|
517 | if (!inserter.second)
|
---|
518 | inserter.first->second = _newvalue;
|
---|
519 | }
|
---|
520 |
|
---|
521 | /** Copies a given trajectory step \a src onto another \a dest
|
---|
522 | * \param dest index of destination step
|
---|
523 | * \param src index of source step
|
---|
524 | */
|
---|
525 | void AtomInfo::CopyStepOnStep(const unsigned int dest, const unsigned int src)
|
---|
526 | {
|
---|
527 | if (dest == src) // self assignment check
|
---|
528 | return;
|
---|
529 |
|
---|
530 | if (WorldTime::getTime() == dest){
|
---|
531 | NOTIFY(AtomObservable::PositionChanged);
|
---|
532 | NOTIFY(AtomObservable::VelocityChanged);
|
---|
533 | NOTIFY(AtomObservable::ForceChanged);
|
---|
534 | }
|
---|
535 |
|
---|
536 | VectorTrajectory_t::iterator positer = AtomicPosition.find(src);
|
---|
537 | ASSERT( positer != AtomicPosition.end(),
|
---|
538 | "AtomInfo::CopyStepOnStep() - step "
|
---|
539 | +toString(src)+" to copy from not present in AtomicPosition.");
|
---|
540 | assignTrajectoryElement(AtomicPosition, dest, positer->second);
|
---|
541 | VectorTrajectory_t::iterator veliter = AtomicVelocity.find(src);
|
---|
542 | if (veliter != AtomicVelocity.end())
|
---|
543 | assignTrajectoryElement(AtomicVelocity, dest, veliter->second);
|
---|
544 | VectorTrajectory_t::iterator forceiter = AtomicForce.find(src);
|
---|
545 | if (forceiter != AtomicForce.end())
|
---|
546 | assignTrajectoryElement(AtomicForce, dest, forceiter->second);
|
---|
547 | };
|
---|
548 |
|
---|
549 | /** Performs a velocity verlet update of the position at \a NextStep from \a LastStep information only.
|
---|
550 | *
|
---|
551 | * We calculate \f$x(t + \delta t) = x(t) + v(t)* \delta t + .5 * \delta t * \delta t * F(t)/m \f$.
|
---|
552 | *
|
---|
553 | *
|
---|
554 | * \param NextStep index of sequential step to set
|
---|
555 | * \param Deltat time step width
|
---|
556 | * \param IsAngstroem whether the force's underlying unit of length is angstroem or bohr radii
|
---|
557 | */
|
---|
558 | void AtomInfo::VelocityVerletUpdateX(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem)
|
---|
559 | {
|
---|
560 | const unsigned int LastStep = NextStep == 0 ? 0 : NextStep-1;
|
---|
561 |
|
---|
562 | LOG(2, "INFO: Particle that currently " << *this);
|
---|
563 | LOG(2, "INFO: Integrating position with mass=" << getMass() << " and Deltat="
|
---|
564 | << Deltat << " at NextStep=" << NextStep);
|
---|
565 |
|
---|
566 | // update position
|
---|
567 | {
|
---|
568 | Vector tempVector = getPositionAtStep(LastStep);
|
---|
569 | LOG(4, "INFO: initial position from last step " << setprecision(4) << tempVector);
|
---|
570 | tempVector += Deltat*(getAtomicVelocityAtStep(LastStep)); // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
|
---|
571 | LOG(4, "INFO: position with velocity " << getAtomicVelocityAtStep(LastStep) << " from last step " << tempVector);
|
---|
572 | tempVector += .5*Deltat*Deltat*(getAtomicForceAtStep(LastStep))*(1./getMass()); // F = m * a and s =
|
---|
573 | LOG(4, "INFO: position with force " << getAtomicForceAtStep(LastStep) << " from last step " << tempVector);
|
---|
574 | setPositionAtStep(NextStep, tempVector);
|
---|
575 | LOG(3, "INFO: Position at step " << NextStep << " set to " << tempVector);
|
---|
576 | }
|
---|
577 | };
|
---|
578 |
|
---|
579 | /** Performs a velocity verlet update of the velocity at \a NextStep.
|
---|
580 | *
|
---|
581 | * \note forces at NextStep should have been calculated based on position at NextStep prior
|
---|
582 | * to calling this function.
|
---|
583 | *
|
---|
584 | * We calculate \f$v(t) = v(t - \delta t) + \delta _t * .5 * (F(t - \delta t) + F(t))/m \f$.
|
---|
585 | *
|
---|
586 | * Parameters are according to those in configuration class.
|
---|
587 | * \param NextStep index of sequential step to set
|
---|
588 | * \param Deltat time step width
|
---|
589 | * \param IsAngstroem whether the force's underlying unit of length is angstroem or bohr radii
|
---|
590 | */
|
---|
591 | void AtomInfo::VelocityVerletUpdateU(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem)
|
---|
592 | {
|
---|
593 | const unsigned int LastStep = NextStep == 0 ? 0 : NextStep-1;
|
---|
594 |
|
---|
595 | LOG(2, "INFO: Particle that currently " << *this);
|
---|
596 | LOG(2, "INFO: Integrating velocity with mass=" << getMass() << " and Deltat="
|
---|
597 | << Deltat << " at NextStep=" << NextStep);
|
---|
598 |
|
---|
599 | // Update U
|
---|
600 | {
|
---|
601 | Vector tempVector = getAtomicVelocityAtStep(LastStep);
|
---|
602 | LOG(4, "INFO: initial velocity from last step " << tempVector);
|
---|
603 | tempVector += Deltat * .5*(getAtomicForceAtStep(LastStep)+getAtomicForceAtStep(NextStep))*(1./getMass()); // v = F/m * t
|
---|
604 | LOG(4, "INFO: Velocity with force from last " << getAtomicForceAtStep(LastStep)
|
---|
605 | << " and present " << getAtomicForceAtStep(NextStep) << " step " << tempVector);
|
---|
606 | setAtomicVelocityAtStep(NextStep, tempVector);
|
---|
607 | LOG(3, "INFO: Velocity at step " << NextStep << " set to " << tempVector);
|
---|
608 | }
|
---|
609 | };
|
---|
610 |
|
---|
611 | std::ostream & AtomInfo::operator << (std::ostream &ost) const
|
---|
612 | {
|
---|
613 | return (ost << getPosition());
|
---|
614 | }
|
---|
615 |
|
---|
616 | std::ostream & operator << (std::ostream &ost, const AtomInfo &a)
|
---|
617 | {
|
---|
618 | const size_t terminalstep = a.getTrajectorySize()-1;
|
---|
619 | if (terminalstep) {
|
---|
620 | ost << "starts at "
|
---|
621 | << a.getPositionAtStep(0) << " and ends at "
|
---|
622 | << a.getPositionAtStep(terminalstep)
|
---|
623 | << " at time step " << terminalstep;
|
---|
624 | } else {
|
---|
625 | ost << "is at "
|
---|
626 | << a.getPositionAtStep(0) << " with a single time step only";
|
---|
627 | }
|
---|
628 | return ost;
|
---|
629 | }
|
---|
630 |
|
---|