Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/Atom/atom_atominfo.cpp

    r8cc22f r7e51e1  
    33 * Description: creates and alters molecular systems
    44 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
    5  * Copyright (C)  2014 Frederik Heber. All rights reserved.
    65 *
    76 *
    … …  
    5655  charge(0.)
    5756{
    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 }
     57  AtomicPosition.reserve(1);
     58  AtomicPosition.push_back(zeroVec);
     59  AtomicVelocity.reserve(1);
     60  AtomicVelocity.push_back(zeroVec);
     61  AtomicForce.reserve(1);
     62  AtomicForce.push_back(zeroVec);
     63
     64};
    6265
    6366/** Copy constructor of class AtomInfo.
    … …  
    6568AtomInfo::AtomInfo(const AtomInfo &_atom) :
    6669    AtomicPosition(_atom.AtomicPosition),
    67     AtomicVelocity(_atom.AtomicVelocity),
    68     AtomicForce(_atom.AtomicForce),
    6970    AtomicElement(_atom.AtomicElement),
    7071    FixedIon(_atom.FixedIon),
    7172    charge(_atom.charge)
    7273{
    73 }
     74  AtomicVelocity.reserve(1);
     75  AtomicVelocity.push_back(zeroVec);
     76  AtomicForce.reserve(1);
     77  AtomicForce.push_back(zeroVec);
     78};
    7479
    7580AtomInfo::AtomInfo(const VectorInterface &_v) :
    … …  
    7883    charge(0.)
    7984{
    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) );
     85  if (AtomicPosition.size() < 1)
     86    AtomicPosition.resize(1);
     87  AtomicPosition[0] = _v.getPosition();
     88  AtomicVelocity.reserve(1);
     89  AtomicVelocity.push_back(zeroVec);
     90  AtomicForce.reserve(1);
     91  AtomicForce.push_back(zeroVec);
    8392};
    8493
    … …  
    8998};
    9099
    91 void AtomInfo::AppendTrajectoryStep(const unsigned int _step)
    92 {
    93   if (WorldTime::getTime() == _step)
    94     NOTIFY(TrajectoryChanged);
    95   AtomicPosition.insert( std::make_pair(_step, zeroVec) );
    96   AtomicVelocity.insert( std::make_pair(_step, zeroVec) );
    97   AtomicForce.insert( std::make_pair(_step, zeroVec) );
     100void AtomInfo::AppendTrajectoryStep()
     101{
     102  NOTIFY(TrajectoryChanged);
     103  AtomicPosition.push_back(zeroVec);
     104  AtomicVelocity.push_back(zeroVec);
     105  AtomicForce.push_back(zeroVec);
    98106  LOG(5,"AtomInfo::AppendTrajectoryStep() called, size is ("
    99107      << AtomicPosition.size() << ","
    … …  
    102110}
    103111
    104 void AtomInfo::removeTrajectoryStep(const unsigned int _step)
    105 {
    106   if (WorldTime::getTime() == _step)
    107     NOTIFY(TrajectoryChanged);
    108   AtomicPosition.erase(_step);
    109   AtomicVelocity.erase(_step);
    110   AtomicForce.erase(_step);
     112void AtomInfo::removeTrajectoryStep()
     113{
     114  NOTIFY(TrajectoryChanged);
     115  AtomicPosition.pop_back();
     116  AtomicVelocity.pop_back();
     117  AtomicForce.pop_back();
    111118  LOG(5,"AtomInfo::removeTrajectoryStep() called, size is ("
    112119      << AtomicPosition.size() << ","
    … …  
    134141const double& AtomInfo::operator[](size_t i) const
    135142{
    136   return atStep(i, WorldTime::getTime());
     143  ASSERT(AtomicPosition.size() > WorldTime::getTime(),
     144      "AtomInfo::operator[]() - Access out of range: "
     145      +toString(WorldTime::getTime())
     146      +" not in [0,"+toString(AtomicPosition.size())+").");
     147  return AtomicPosition[WorldTime::getTime()][i];
    137148}
    138149
    139150const double& AtomInfo::at(size_t i) const
    140151{
    141   return atStep(i, WorldTime::getTime());
     152  ASSERT(AtomicPosition.size() > WorldTime::getTime(),
     153      "AtomInfo::at() - Access out of range: "
     154      +toString(WorldTime::getTime())
     155      +" not in [0,"+toString(AtomicPosition.size())+").");
     156  return AtomicPosition[WorldTime::getTime()].at(i);
    142157}
    143158
    144159const double& AtomInfo::atStep(size_t i, unsigned int _step) const
    145160{
    146   ASSERT(!AtomicPosition.empty(),
    147       "AtomInfo::operator[]() - AtomicPosition is empty.");
    148   VectorTrajectory_t::const_iterator iter =
    149       AtomicPosition.lower_bound(_step);
    150   return iter->second[i];
     161  ASSERT(AtomicPosition.size() > _step,
     162      "AtomInfo::atStep() - Access out of range: "
     163      +toString(_step)
     164      +" not in [0,"+toString(AtomicPosition.size())+").");
     165  return AtomicPosition[_step].at(i);
    151166}
    152167
    … …  
    155170  OBSERVE;
    156171  NOTIFY(AtomObservable::PositionChanged);
    157   VectorTrajectory_t::iterator iter = AtomicPosition.find(WorldTime::getTime());
    158   if (iter !=  AtomicPosition.end()) {
    159     iter->second[i] = value;
    160   } else {
    161     Vector newPos;
    162     newPos[i] = value;
    163 #ifndef NDEBUG
    164     std::pair<VectorTrajectory_t::iterator, bool> inserter =
    165 #endif
    166         AtomicPosition.insert( std::make_pair(WorldTime::getTime(), newPos) );
    167     ASSERT( inserter.second,
    168         "AtomInfo::set() - time step "+toString(WorldTime::getTime())
    169         +" present after all?");
    170   }
    171 }
    172 
    173 /** Helps to determine whether the current step really exists or getPosition() has just
    174  * delivered the one closest to it in the past.
    175  *
    176  * \param _step step to check
    177  * \param true - step exists, false - step does not exist, getPosition() delivers closest
    178  */
    179 bool AtomInfo::isStepPresent(const unsigned int _step) const
    180 {
    181   VectorTrajectory_t::const_iterator iter =
    182       AtomicPosition.find(_step);
    183   return iter != AtomicPosition.end();
     172  ASSERT(AtomicPosition.size() > WorldTime::getTime(),
     173      "AtomInfo::set() - Access out of range: "
     174      +toString(WorldTime::getTime())
     175      +" not in [0,"+toString(AtomicPosition.size())+").");
     176  AtomicPosition[WorldTime::getTime()].at(i) = value;
    184177}
    185178
    186179const Vector& AtomInfo::getPosition() const
    187180{
    188   return getPositionAtStep(WorldTime::getTime());
     181  ASSERT(AtomicPosition.size() > WorldTime::getTime(),
     182      "AtomInfo::getPosition() - Access out of range: "
     183      +toString(WorldTime::getTime())
     184      +" not in [0,"+toString(AtomicPosition.size())+").");
     185  return AtomicPosition[WorldTime::getTime()];
    189186}
    190187
    191188const Vector& AtomInfo::getPositionAtStep(const unsigned int _step) const
    192189{
    193   ASSERT(!AtomicPosition.empty(),
    194       "AtomInfo::operator[]() - AtomicPosition is empty.");
    195   VectorTrajectory_t::const_iterator iter =
    196       AtomicPosition.lower_bound(_step);
    197   return iter->second;
     190  ASSERT(_step < AtomicPosition.size(),
     191      "AtomInfo::getPositionAtStep() - Access out of range: "
     192      +toString(_step)
     193      +" not in [0,"+toString(AtomicPosition.size())+").");
     194  return AtomicPosition[_step];
    198195}
    199196
    … …  
    210207{
    211208  const element *elem = World::getInstance().getPeriode()->FindElement(Z);
    212   setType(elem);
    213 }
     209  if (elem != NULL) {
     210    OBSERVE;
     211    NOTIFY(AtomObservable::ElementChanged);
     212    AtomicElement = Z;
     213  }
     214}
     215
     216//Vector& AtomInfo::getAtomicVelocity()
     217//{
     218//  return AtomicVelocity[0];
     219//}
     220
     221//Vector& AtomInfo::getAtomicVelocity(const int _step)
     222//{
     223//  ASSERT(_step < AtomicVelocity.size(),
     224//      "AtomInfo::getAtomicVelocity() - Access out of range.");
     225//  return AtomicVelocity[_step];
     226//}
    214227
    215228const Vector& AtomInfo::getAtomicVelocity() const
    216229{
    217   return getAtomicVelocityAtStep(WorldTime::getTime());
     230  ASSERT(AtomicVelocity.size() > 0,
     231      "AtomInfo::getAtomicVelocity() - Access out of range: "
     232      +toString(WorldTime::getTime())
     233      +" not in [0,"+toString(AtomicPosition.size())+").");
     234  return AtomicVelocity[WorldTime::getTime()];
    218235}
    219236
    220237const Vector& AtomInfo::getAtomicVelocityAtStep(const unsigned int _step) const
    221238{
    222   ASSERT(!AtomicVelocity.empty(),
    223       "AtomInfo::operator[]() - AtomicVelocity is empty.");
    224   VectorTrajectory_t::const_iterator iter =
    225       AtomicVelocity.lower_bound(_step);
    226   // special, we only interpolate between present time steps not into the future
    227   if (_step > AtomicVelocity.begin()->first)
    228     return zeroVec;
    229   else
    230     return iter->second;
     239  ASSERT(_step < AtomicVelocity.size(),
     240      "AtomInfo::getAtomicVelocity() - Access out of range: "
     241      +toString(_step)
     242      +" not in [0,"+toString(AtomicPosition.size())+").");
     243  return AtomicVelocity[_step];
    231244}
    232245
    233246void AtomInfo::setAtomicVelocity(const Vector &_newvelocity)
    234247{
    235   setAtomicVelocityAtStep(WorldTime::getTime(), _newvelocity);
     248  OBSERVE;
     249  NOTIFY(AtomObservable::VelocityChanged);
     250  ASSERT(WorldTime::getTime() < AtomicVelocity.size(),
     251      "AtomInfo::setAtomicVelocity() - Access out of range: "
     252      +toString(WorldTime::getTime())
     253      +" not in [0,"+toString(AtomicPosition.size())+").");
     254  AtomicVelocity[WorldTime::getTime()] = _newvelocity;
    236255}
    237256
    … …  
    239258{
    240259  OBSERVE;
    241   VectorTrajectory_t::iterator iter = AtomicVelocity.find(_step);
    242   if (iter !=  AtomicVelocity.end()) {
    243     iter->second = _newvelocity;
    244   } else {
    245 #ifndef NDEBUG
    246     std::pair<VectorTrajectory_t::iterator, bool> inserter =
    247 #endif
    248         AtomicVelocity.insert( std::make_pair(_step, _newvelocity) );
    249     ASSERT( inserter.second,
    250         "AtomInfo::set() - time step "+toString(_step)
    251         +" present after all?");
    252   }
    253260  if (WorldTime::getTime() == _step)
    254261    NOTIFY(AtomObservable::VelocityChanged);
     262  const unsigned int size = AtomicVelocity.size();
     263  ASSERT(_step <= size,
     264      "AtomInfo::setAtomicVelocityAtStep() - Access out of range: "
     265      +toString(_step)
     266      +" not in [0,"+toString(size)+"].");
     267  if(_step < size) {
     268    AtomicVelocity[_step] = _newvelocity;
     269  } else if (_step == size) {
     270    UpdateSteps();
     271    AtomicVelocity[_step] = _newvelocity;
     272  }
    255273}
    256274
    257275const Vector& AtomInfo::getAtomicForce() const
    258276{
    259   return getAtomicForceAtStep(WorldTime::getTime());
     277  ASSERT(WorldTime::getTime() < AtomicForce.size(),
     278      "AtomInfo::getAtomicForce() - Access out of range: "
     279      +toString(WorldTime::getTime())
     280      +" not in [0,"+toString(AtomicPosition.size())+").");
     281  return AtomicForce[WorldTime::getTime()];
    260282}
    261283
    262284const Vector& AtomInfo::getAtomicForceAtStep(const unsigned int _step) const
    263285{
    264   ASSERT(!AtomicForce.empty(),
    265       "AtomInfo::operator[]() - AtomicForce is empty.");
    266   VectorTrajectory_t::const_iterator iter =
    267       AtomicForce.lower_bound(_step);
    268   // special, we only interpolate between present time steps not into the future
    269   if (_step > AtomicForce.begin()->first)
    270     return zeroVec;
    271   else
    272     return iter->second;
     286  ASSERT(_step < AtomicForce.size(),
     287      "AtomInfo::getAtomicForce() - Access out of range: "
     288      +toString(_step)
     289      +" not in [0,"+toString(AtomicPosition.size())+").");
     290  return AtomicForce[_step];
    273291}
    274292
    275293void AtomInfo::setAtomicForce(const Vector &_newforce)
    276294{
    277   setAtomicForceAtStep(WorldTime::getTime(), _newforce);
     295  OBSERVE;
     296  NOTIFY(AtomObservable::ForceChanged);
     297  ASSERT(WorldTime::getTime() < AtomicForce.size(),
     298      "AtomInfo::setAtomicForce() - Access out of range: "
     299      +toString(WorldTime::getTime())
     300      +" not in [0,"+toString(AtomicPosition.size())+").");
     301  AtomicForce[WorldTime::getTime()] = _newforce;
    278302}
    279303
    … …  
    281305{
    282306  OBSERVE;
    283   VectorTrajectory_t::iterator iter = AtomicForce.find(_step);
    284   if (iter !=  AtomicForce.end()) {
    285     iter->second = _newforce;
    286   } else {
    287 #ifndef NDEBUG
    288     std::pair<VectorTrajectory_t::iterator, bool> inserter =
    289 #endif
    290         AtomicForce.insert( std::make_pair(_step, _newforce) );
    291     ASSERT( inserter.second,
    292         "AtomInfo::set() - time step "+toString(_step)
    293         +" present after all?");
    294   }
    295307  if (WorldTime::getTime() == _step)
    296308    NOTIFY(AtomObservable::ForceChanged);
     309  const unsigned int size = AtomicForce.size();
     310  ASSERT(_step <= size,
     311      "AtomInfo::setAtomicForce() - Access out of range: "
     312      +toString(_step)
     313      +" not in [0,"+toString(AtomicPosition.size())+"].");
     314  if(_step < size) {
     315    AtomicForce[_step] = _newforce;
     316  } else if (_step == size) {
     317    UpdateSteps();
     318    AtomicForce[_step] = _newforce;
     319  }
    297320}
    298321
    … …  
    311334void AtomInfo::setPosition(const Vector& _vector)
    312335{
    313   setPositionAtStep(WorldTime::getTime(), _vector);
     336  OBSERVE;
     337  NOTIFY(AtomObservable::PositionChanged);
     338  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     339      "AtomInfo::setPosition() - Access out of range: "
     340      +toString(WorldTime::getTime())
     341      +" not in [0,"+toString(AtomicPosition.size())+").");
     342  AtomicPosition[WorldTime::getTime()] = _vector;
     343  //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
    314344}
    315345
    … …  
    317347{
    318348  OBSERVE;
    319   VectorTrajectory_t::iterator iter = AtomicPosition.find(_step);
    320   if (iter !=  AtomicPosition.end()) {
    321     iter->second = _vector;
    322   } else {
    323 #ifndef NDEBUG
    324     std::pair<VectorTrajectory_t::iterator, bool> inserter =
    325 #endif
    326         AtomicPosition.insert( std::make_pair(_step, _vector) );
    327     ASSERT( inserter.second,
    328         "AtomInfo::set() - time step "+toString(_step)
    329         +" present after all?");
    330   }
    331349  if (WorldTime::getTime() == _step)
    332350    NOTIFY(AtomObservable::PositionChanged);
     351  const unsigned int size = AtomicPosition.size();
     352  ASSERT(_step <= size,
     353      "AtomInfo::setPosition() - Access out of range: "
     354      +toString(_step)
     355      +" not in [0,"+toString(size)+"].");
     356  if(_step < size) {
     357    AtomicPosition[_step] = _vector;
     358  } else if (_step == size) {
     359    UpdateSteps();
     360    AtomicPosition[_step] = _vector;
     361  }
     362  //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
    333363}
    334364
    335365const VectorInterface& AtomInfo::operator+=(const Vector& b)
    336366{
    337   setPosition(getPosition()+b);
     367  OBSERVE;
     368  NOTIFY(AtomObservable::PositionChanged);
     369  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     370      "AtomInfo::operator+=() - Access out of range: "
     371      +toString(WorldTime::getTime())
     372      +" not in [0,"+toString(AtomicPosition.size())+").");
     373  AtomicPosition[WorldTime::getTime()] += b;
    338374  return *this;
    339375}
    … …  
    341377const VectorInterface& AtomInfo::operator-=(const Vector& b)
    342378{
    343   setPosition(getPosition()-b);
     379  OBSERVE;
     380  NOTIFY(AtomObservable::PositionChanged);
     381  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     382      "AtomInfo::operator-=() - Access out of range: "
     383      +toString(WorldTime::getTime())
     384      +" not in [0,"+toString(AtomicPosition.size())+").");
     385  AtomicPosition[WorldTime::getTime()] -= b;
    344386  return *this;
    345387}
    … …  
    347389Vector const AtomInfo::operator+(const Vector& b) const
    348390{
    349   Vector a(getPosition());
     391  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     392      "AtomInfo::operator+() - Access out of range: "
     393      +toString(WorldTime::getTime())
     394      +" not in [0,"+toString(AtomicPosition.size())+").");
     395  Vector a(AtomicPosition[WorldTime::getTime()]);
    350396  a += b;
    351397  return a;
    … …  
    354400Vector const AtomInfo::operator-(const Vector& b) const
    355401{
    356   Vector a(getPosition());
     402  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     403      "AtomInfo::operator-() - Access out of range: "
     404      +toString(WorldTime::getTime())
     405      +" not in [0,"+toString(AtomicPosition.size())+").");
     406  Vector a(AtomicPosition[WorldTime::getTime()]);
    357407  a -= b;
    358408  return a;
    … …  
    361411double AtomInfo::distance(const Vector &point) const
    362412{
    363   return getPosition().distance(point);
     413  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     414      "AtomInfo::distance() - Access out of range: "
     415      +toString(WorldTime::getTime())
     416      +" not in [0,"+toString(AtomicPosition.size())+").");
     417  return AtomicPosition[WorldTime::getTime()].distance(point);
    364418}
    365419
    366420double AtomInfo::DistanceSquared(const Vector &y) const
    367421{
    368   return getPosition().DistanceSquared(y);
     422  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     423      "AtomInfo::DistanceSquared() - Access out of range: "
     424      +toString(WorldTime::getTime())
     425      +" not in [0,"+toString(AtomicPosition.size())+").");
     426  return AtomicPosition[WorldTime::getTime()].DistanceSquared(y);
    369427}
    370428
    371429double AtomInfo::distance(const VectorInterface &_atom) const
    372430{
    373   return _atom.distance(getPosition());
     431  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     432      "AtomInfo::distance() - Access out of range: "
     433      +toString(WorldTime::getTime())
     434      +" not in [0,"+toString(AtomicPosition.size())+").");
     435  return _atom.distance(AtomicPosition[WorldTime::getTime()]);
    374436}
    375437
    376438double AtomInfo::DistanceSquared(const VectorInterface &_atom) const
    377439{
    378   return _atom.DistanceSquared(getPosition());
     440  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     441      "AtomInfo::DistanceSquared() - Access out of range: "
     442      +toString(WorldTime::getTime())
     443      +" not in [0,"+toString(AtomicPosition.size())+").");
     444  return _atom.DistanceSquared(AtomicPosition[WorldTime::getTime()]);
    379445}
    380446
    381447VectorInterface &AtomInfo::operator=(const Vector& _vector)
    382448{
    383   setPosition(_vector);
     449  OBSERVE;
     450  NOTIFY(AtomObservable::PositionChanged);
     451  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     452      "AtomInfo::operator=() - Access out of range: "
     453      +toString(WorldTime::getTime())
     454      +" not in [0,"+toString(AtomicPosition.size())+").");
     455  AtomicPosition[WorldTime::getTime()] = _vector;
    384456  return *this;
    385457}
    … …  
    387459void AtomInfo::ScaleAll(const double *factor)
    388460{
    389   Vector temp(getPosition());
    390   temp.ScaleAll(factor);
    391   setPosition(temp);
     461  OBSERVE;
     462  NOTIFY(AtomObservable::PositionChanged);
     463  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     464      "AtomInfo::ScaleAll() - Access out of range: "
     465      +toString(WorldTime::getTime())
     466      +" not in [0,"+toString(AtomicPosition.size())+").");
     467  AtomicPosition[WorldTime::getTime()].ScaleAll(factor);
    392468}
    393469
    394470void AtomInfo::ScaleAll(const Vector &factor)
    395471{
    396   Vector temp(getPosition());
    397   temp.ScaleAll(factor);
    398   setPosition(temp);
     472  OBSERVE;
     473  NOTIFY(AtomObservable::PositionChanged);
     474  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     475      "AtomInfo::ScaleAll() - Access out of range: "
     476      +toString(WorldTime::getTime())
     477      +" not in [0,"+toString(AtomicPosition.size())+").");
     478  AtomicPosition[WorldTime::getTime()].ScaleAll(factor);
    399479}
    400480
    401481void AtomInfo::Scale(const double factor)
    402482{
    403   Vector temp(getPosition());
    404   temp.Scale(factor);
    405   setPosition(temp);
     483  OBSERVE;
     484  NOTIFY(AtomObservable::PositionChanged);
     485  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     486      "AtomInfo::Scale() - Access out of range: "
     487      +toString(WorldTime::getTime())
     488      +" not in [0,"+toString(AtomicPosition.size())+").");
     489  AtomicPosition[WorldTime::getTime()].Scale(factor);
    406490}
    407491
    408492void AtomInfo::Zero()
    409493{
    410   setPosition(zeroVec);
     494  OBSERVE;
     495  NOTIFY(AtomObservable::PositionChanged);
     496  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     497      "AtomInfo::Zero() - Access out of range: "
     498      +toString(WorldTime::getTime())
     499      +" not in [0,"+toString(AtomicPosition.size())+").");
     500  AtomicPosition[WorldTime::getTime()].Zero();
    411501}
    412502
    413503void AtomInfo::One(const double one)
    414504{
    415   setPosition(Vector(one,one,one));
     505  OBSERVE;
     506  NOTIFY(AtomObservable::PositionChanged);
     507  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     508      "AtomInfo::One() - Access out of range: "
     509      +toString(WorldTime::getTime())
     510      +" not in [0,"+toString(AtomicPosition.size())+").");
     511  AtomicPosition[WorldTime::getTime()].One(one);
    416512}
    417513
    418514void AtomInfo::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
    419515{
    420   Vector newPos;
    421   newPos.LinearCombinationOfVectors(x1,x2,x3,factors);
    422   setPosition(newPos);
     516  OBSERVE;
     517  NOTIFY(AtomObservable::PositionChanged);
     518  ASSERT(WorldTime::getTime() < AtomicPosition.size(),
     519      "AtomInfo::LinearCombinationOfVectors() - Access out of range: "
     520      +toString(WorldTime::getTime())
     521      +" not in [0,"+toString(AtomicPosition.size())+").");
     522  AtomicPosition[WorldTime::getTime()].LinearCombinationOfVectors(x1,x2,x3,factors);
    423523}
    424524
    … …  
    428528double AtomInfo::getKineticEnergy(const unsigned int _step) const
    429529{
    430   return getMass() * getAtomicVelocityAtStep(_step).NormSquared();
     530  ASSERT(_step < AtomicPosition.size(),
     531      "AtomInfo::getKineticEnergy() - Access out of range: "
     532      +toString(WorldTime::getTime())
     533      +" not in [0,"+toString(AtomicPosition.size())+").");
     534  return getMass() * AtomicVelocity[_step].NormSquared();
    431535}
    432536
    433537Vector AtomInfo::getMomentum(const unsigned int _step) const
    434538{
    435   return getMass() * getAtomicVelocityAtStep(_step);
    436 }
    437 
    438 /** Decrease the trajectory if given \a MaxSteps is smaller.
    439  * Does nothing if \a MaxSteps is larger than current size.
    440  *
     539  ASSERT(_step < AtomicPosition.size(),
     540      "AtomInfo::getMomentum() - Access out of range: "
     541      +toString(WorldTime::getTime())
     542      +" not in [0,"+toString(AtomicPosition.size())+").");
     543  return getMass()*AtomicVelocity[_step];
     544}
     545
     546/** Extends the trajectory STL vector to the new size.
     547 * Does nothing if \a MaxSteps is smaller than current size.
    441548 * \param MaxSteps
    442549 */
    443550void AtomInfo::ResizeTrajectory(size_t MaxSteps)
    444551{
    445   // mind the reverse ordering due to std::greater, latest time steps are at beginning
    446   VectorTrajectory_t::iterator positer = AtomicPosition.lower_bound(MaxSteps);
    447   if (positer != AtomicPosition.begin()) {
    448     if (positer->first == MaxSteps)
    449       --positer;
    450     AtomicPosition.erase(AtomicPosition.begin(), positer);
    451   }
    452   VectorTrajectory_t::iterator veliter = AtomicVelocity.lower_bound(MaxSteps);
    453   if (veliter != AtomicVelocity.begin()) {
    454     if (veliter->first == MaxSteps)
    455       --veliter;
    456     AtomicVelocity.erase(AtomicVelocity.begin(), veliter);
    457   }
    458   VectorTrajectory_t::iterator forceiter = AtomicForce.lower_bound(MaxSteps);
    459   if (forceiter != AtomicForce.begin()) {
    460     if (forceiter->first == MaxSteps)
    461       --forceiter;
    462     AtomicForce.erase(AtomicForce.begin(), forceiter);
    463   }
     552  for (;AtomicPosition.size() <= (unsigned int)(MaxSteps);)
     553    UpdateSteps();
    464554}
    465555
    466556size_t AtomInfo::getTrajectorySize() const
    467557{
    468   // mind greater comp for map here: first element is latest in time steps!
    469   return AtomicPosition.begin()->first+1;
     558  return AtomicPosition.size();
    470559}
    471560
    … …  
    473562{
    474563  return getType()->getMass();
    475 }
    476 
    477 /** Helper function to either insert or assign, depending on the element being
    478  * present already.
    479  *
    480  * \param _trajectory vector of Vectors to assign
    481  * \param dest step to insert/assign to
    482  * \param _newvalue new Vector value
    483  */
    484 void assignTrajectoryElement(
    485     std::map<unsigned int, Vector, std::greater<unsigned int> > &_trajectory,
    486     const unsigned int dest,
    487     const Vector &_newvalue)
    488 {
    489   std::pair<std::map<unsigned int, Vector, std::greater<unsigned int> >::iterator, bool> inserter =
    490       _trajectory.insert( std::make_pair(dest, _newvalue) );
    491   if (!inserter.second)
    492     inserter.first->second = _newvalue;
    493564}
    494565
    … …  
    508579  }
    509580
    510   VectorTrajectory_t::iterator positer = AtomicPosition.find(src);
    511   ASSERT( positer != AtomicPosition.end(),
    512       "AtomInfo::CopyStepOnStep() - step "
    513       +toString(src)+" to copy from not present in AtomicPosition.");
    514   VectorTrajectory_t::iterator veliter = AtomicVelocity.find(src);
    515   ASSERT( veliter != AtomicVelocity.end(),
    516       "AtomInfo::CopyStepOnStep() - step "
    517       +toString(src)+" to copy from not present in AtomicVelocity.");
    518   VectorTrajectory_t::iterator forceiter = AtomicForce.find(src);
    519   ASSERT( forceiter != AtomicForce.end(),
    520       "AtomInfo::CopyStepOnStep() - step "
    521       +toString(src)+" to copy from not present in AtomicForce.");
    522   assignTrajectoryElement(AtomicPosition, dest, positer->second);
    523   assignTrajectoryElement(AtomicVelocity, dest, veliter->second);
    524   assignTrajectoryElement(AtomicForce, dest, forceiter->second);
     581  ASSERT(dest < AtomicPosition.size(),
     582      "AtomInfo::CopyStepOnStep() - destination outside of current trajectory array: "
     583      +toString(dest)
     584      +" not in [0,"+toString(AtomicPosition.size())+").");
     585  ASSERT(src < AtomicPosition.size(),
     586      "AtomInfo::CopyStepOnStep() - source outside of current trajectory array: "
     587      +toString(src)
     588      +" not in [0,"+toString(AtomicPosition.size())+").");
     589  for (int n=NDIM;n--;) {
     590    AtomicPosition.at(dest)[n] = AtomicPosition.at(src)[n];
     591    AtomicVelocity.at(dest)[n] = AtomicVelocity.at(src)[n];
     592    AtomicForce.at(dest)[n] = AtomicForce.at(src)[n];
     593  }
    525594};
    526595
    … …  
    587656};
    588657
     658//const AtomInfo& operator*=(AtomInfo& a, const double m)
     659//{
     660//  a.Scale(m);
     661//  return a;
     662//}
     663//
     664//AtomInfo const operator*(const AtomInfo& a, const double m)
     665//{
     666//  AtomInfo copy(a);
     667//  copy *= m;
     668//  return copy;
     669//}
     670//
     671//AtomInfo const operator*(const double m, const AtomInfo& a)
     672//{
     673//  AtomInfo copy(a);
     674//  copy *= m;
     675//  return copy;
     676//}
     677
    589678std::ostream & AtomInfo::operator << (std::ostream &ost) const
    590679{
Note: See TracChangeset for help on using the changeset viewer.