Changeset 056e70


Ignore:
Timestamp:
Feb 24, 2011, 5:56:03 PM (14 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
6b020f
Parents:
6625c3
Message:

Suffixed getters and setters for AtomInfo trajecories with AtStep.

  • AtomInfo:: getters are now const members
  • added lots of ASSERTS to getters with time step to assure no out-of-bounds access.
Location:
src
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • src/AtomSet.hpp

    r6625c3 r056e70  
    127127inline void AtomSetMixin<Set>::addVelocityAtStep(const Vector velocity, unsigned int step){
    128128  BOOST_FOREACH(AtomInfo *atom,*this){
    129     atom->getAtomicVelocity(step) += velocity;
     129    atom->setAtomicVelocityAtStep(step, atom->getAtomicVelocityAtStep(step)+velocity);
    130130  }
    131131}

  • src/Parser/TremoloParser.cpp

    r6625c3 r056e70  
    313313  ConvertTo<double> toDouble;
    314314  ConvertTo<int> toInt;
     315  Vector tempVector;
    315316
    316317  // setup tokenizer, splitting up white-spaced entries
    … …  
    327328    currentField = knownKeys[keyName];
    328329    const string word = *tok_iter;
    329     //DoLog(1) && (Log() << Verbose(1) << "INFO: Parsing key " << keyName << " with remaining data " << word << std::endl);
     330    DoLog(4) && (Log() << Verbose(4) << "INFO: Parsing key " << keyName << " with remaining data " << word << std::endl);
    330331    switch (currentField) {
    331332      case TremoloKey::x :
    … …  
    333334        for (int i=0;i<NDIM;i++) {
    334335          ASSERT(tok_iter != tokens.end(), "TremoloParser::readAtomDataLine() - no value for x["+toString(i)+"]!");
    335           //DoLog(1) && (Log() << Verbose(1) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
     336          DoLog(4) && (Log() << Verbose(4) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
    336337          newAtom->set(i, toDouble(*tok_iter));
    337338          tok_iter++;
    … …  
    342343        for (int i=0;i<NDIM;i++) {
    343344          ASSERT(tok_iter != tokens.end(), "TremoloParser::readAtomDataLine() - no value for u["+toString(i)+"]!");
    344           //DoLog(1) && (Log() << Verbose(1) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
    345           newAtom->getAtomicVelocity()[i] = toDouble(*tok_iter);
     345          DoLog(4) && (Log() << Verbose(4) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
     346          tempVector[i] = toDouble(*tok_iter);
    346347          tok_iter++;
    347348        }
     349        newAtom->setAtomicVelocity(tempVector);
    348350        break;
    349351      case TremoloKey::Type :
    350352      {
    351353        ASSERT(tok_iter != tokens.end(), "TremoloParser::readAtomDataLine() - no value for "+keyName+"!");
    352         DoLog(1) && (Log() << Verbose(1) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
     354        DoLog(4) && (Log() << Verbose(4) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
    353355        std::string element(knownTypes[(*tok_iter)]);
    354356        // put type name into container for later use
    355357        atomInfo->set(currentField, *tok_iter);
    356         DoLog(1) && (Log() << Verbose(1) << "INFO: Parsing element " << (*tok_iter) << " as " << element << " according to KnownTypes." << std::endl);
     358        DoLog(4) && (Log() << Verbose(4) << "INFO: Parsing element " << (*tok_iter) << " as " << element << " according to KnownTypes." << std::endl);
    357359        tok_iter++;
    358360        newAtom->setType(World::getInstance().getPeriode()->FindElement(element));
    … …  
    362364      case TremoloKey::Id :
    363365        ASSERT(tok_iter != tokens.end(), "TremoloParser::readAtomDataLine() - no value for "+keyName+"!");
    364         //DoLog(1) && (Log() << Verbose(1) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
     366        DoLog(4) && (Log() << Verbose(4) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
    365367        atomIdMap[toInt(*tok_iter)] = newAtom->getId();
    366368        tok_iter++;
    … …  
    369371        for (int i=0;i<atoi(it->substr(it->find("=") + 1, 1).c_str());i++) {
    370372          ASSERT(tok_iter != tokens.end(), "TremoloParser::readAtomDataLine() - no value for "+keyName+"!");
    371           //DoLog(1) && (Log() << Verbose(1) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
     373          DoLog(4) && (Log() << Verbose(4) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
    372374          lineStream << *tok_iter << "\t";
    373375          tok_iter++;
    … …  
    378380      default :
    379381        ASSERT(tok_iter != tokens.end(), "TremoloParser::readAtomDataLine() - no value for "+keyName+"!");
    380         //DoLog(1) && (Log() << Verbose(1) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
     382        DoLog(4) && (Log() << Verbose(4) << "INFO: Parsing key " << keyName << " with next token " << *tok_iter << std::endl);
    381383        atomInfo->set(currentField, *tok_iter);
    382384        tok_iter++;
    … …  
    403405    // 0 is used to fill empty neighbor positions in the tremolo file.
    404406    if (neighborId > 0) {
    405 //      DoLog(1) && (Log() << Verbose(1)
    406 //          << "Atom with global id " << atomId
    407 //          << " has neighbour with serial " << neighborId
    408 //          << std::endl);
     407      DoLog(4) && (Log() << Verbose(4)
     408          << "Atom with global id " << atomId
     409          << " has neighbour with serial " << neighborId
     410          << std::endl);
    409411      additionalAtomData[atomId].neighbors.push_back(neighborId);
    410412    }

  • src/Thermostats/Berendsen.cpp

    r6625c3 r056e70  
    6363  double ScaleTempFactor = getContainer().TargetTemp/ActualTemp;
    6464  for(ForwardIterator iter=begin;iter!=end;++iter){
    65     Vector &U = (*iter)->getAtomicVelocity(step);
     65    Vector U = (*iter)->getAtomicVelocityAtStep(step);
    6666    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
    6767      U *= sqrt(1+(World::getInstance().getConfig()->Deltat/TempFrequency)*(ScaleTempFactor-1));
    6868      ekin += 0.5*(*iter)->getType()->getMass() * U.NormSquared();
    6969    }
     70    (*iter)->setAtomicVelocityAtStep(step, U);
    7071  }
    7172  return ekin;

  • src/Thermostats/GaussianThermostat.cpp

    r6625c3 r056e70  
    7070  double ekin =0;
    7171  for(ForwardIterator iter=begin;iter!=end;++iter){
    72     Vector &U = (*iter)->getAtomicVelocity(step);
     72    Vector U = (*iter)->getAtomicVelocityAtStep(step);
    7373    if ((*iter)->getFixedIon() == 0) {// even FixedIon moves, only not by other's forces
    7474      U += World::getInstance().getConfig()->Deltat * G_over_E * U;
    7575      ekin += (*iter)->getType()->getMass() * U.NormSquared();
    7676    }
     77    (*iter)->setAtomicVelocityAtStep(step, U);
    7778  }
    7879  return ekin;
    … …  
    8485  G=0;
    8586  for(ForwardIterator iter=begin;iter!=end;++iter){
    86     const Vector &U = (*iter)->getAtomicVelocity(step);
    87     const Vector &F = (*iter)->getAtomicForce(step);
     87    const Vector &U = (*iter)->getAtomicVelocityAtStep(step);
     88    const Vector &F = (*iter)->getAtomicForceAtStep(step);
    8889    if ((*iter)->getFixedIon() == 0){ // even FixedIon moves, only not by other's forces
    8990      G += U.ScalarProduct(F);

  • src/Thermostats/Langevin.cpp

    r6625c3 r056e70  
    8989    rng_distribution = new boost::normal_distribution<>(0,sigma);
    9090    boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > rng(*rng_engine, *rng_distribution);
    91     Vector &U = (*iter)->getAtomicVelocity(step);
     91    Vector U = (*iter)->getAtomicVelocityAtStep(step);
    9292    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
    9393      // throw a dice to determine whether it gets hit by a heat bath particle
    … …  
    102102      ekin += 0.5*(*iter)->getType()->getMass() * U.NormSquared();
    103103    }
     104    (*iter)->setAtomicVelocityAtStep(step, U);
    104105    delete rng_distribution;
    105106    rng_distribution = NULL;

  • src/Thermostats/NoseHoover.cpp

    r6625c3 r056e70  
    6868  double ekin =0;
    6969  for(ForwardIterator iter=begin;iter!=end;++iter){
    70     Vector &U = (*iter)->getAtomicVelocity(step);
     70    Vector U = (*iter)->getAtomicVelocityAtStep(step);
    7171    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
    7272      U += World::getInstance().getConfig()->Deltat/(*iter)->getType()->getMass() * (alpha * (U * (*iter)->getType()->getMass()));
    7373      ekin += (0.5*(*iter)->getType()->getMass()) * U.NormSquared();
    7474    }
     75    (*iter)->setAtomicVelocityAtStep(step, U);
    7576  }
    7677  return ekin;
    … …  
    8283  count=0;
    8384  for(ForwardIterator iter = begin;iter!=end;++iter){
    84     const Vector &U = (*iter)->getAtomicVelocity(step);
     85    const Vector &U = (*iter)->getAtomicVelocityAtStep(step);
    8586    if ((*iter)->getFixedIon() == 0) { // even FixedIon moves, only not by other's forces
    8687      delta_alpha += U.NormSquared()*(*iter)->getType()->getMass();

  • src/Thermostats/Woodcock.cpp

    r6625c3 r056e70  
    6767    double ekin;
    6868    for (ForwardIterator iter = begin; iter!=end;++iter){
    69       Vector &U = (*iter)->getAtomicVelocity(step);
     69      Vector U = (*iter)->getAtomicVelocityAtStep(step);
    7070      if ((*iter)->getFixedIon() == 0){ // even FixedIon moves, only not by other's forces
    7171        U *= ScaleTempFactor;
    7272        ekin += 0.5*(*iter)->getType()->getMass() * U.NormSquared();
    7373      }
     74      (*iter)->setAtomicVelocityAtStep(step, U);
    7475    }
    7576  }

  • src/atom.cpp

    r6625c3 r056e70  
    218218    ASSERT(elementLookup.there.find(elemental)!=elementLookup.there.end(),"Type of this atom was not in the formula upon enumeration");
    219219    *out << "Ion_Type" << elementLookup.there.find(elemental)->second << "_" << AtomNo[getType()->getAtomicNumber()] << "\t"  << fixed << setprecision(9) << showpoint;
    220     *out << getPosition(step)[0] << "\t" << getPosition(step)[1] << "\t" << getPosition(step)[2];
     220    *out << getPositionAtStep(step)[0] << "\t" << getPositionAtStep(step)[1] << "\t" << getPositionAtStep(step)[2];
    221221    *out << "\t" << (int)(getFixedIon());
    222     if (getAtomicVelocity(step).Norm() > MYEPSILON)
    223       *out << "\t" << scientific << setprecision(6) << getAtomicVelocity(step)[0] << "\t" << getAtomicVelocity(step)[1] << "\t" << getAtomicVelocity(step)[2] << "\t";
    224     if (getAtomicForce(step).Norm() > MYEPSILON)
    225       *out << "\t" << scientific << setprecision(6) << getAtomicForce(step)[0] << "\t" << getAtomicForce(step)[1] << "\t" << getAtomicForce(step)[2] << "\t";
     222    if (getAtomicVelocityAtStep(step).Norm() > MYEPSILON)
     223      *out << "\t" << scientific << setprecision(6) << getAtomicVelocityAtStep(step)[0] << "\t" << getAtomicVelocityAtStep(step)[1] << "\t" << getAtomicVelocityAtStep(step)[2] << "\t";
     224    if (getAtomicForceAtStep(step).Norm() > MYEPSILON)
     225      *out << "\t" << scientific << setprecision(6) << getAtomicForceAtStep(step)[0] << "\t" << getAtomicForceAtStep(step)[1] << "\t" << getAtomicForceAtStep(step)[2] << "\t";
    226226    *out << "\t# Number in molecule " << nr << endl;
    227227    return true;
    … …  
    239239  if (out != NULL) {
    240240    *out << getType()->getSymbol() << "\t";
    241     *out << getPosition(step)[0] << "\t";
    242     *out << getPosition(step)[1] << "\t";
    243     *out << getPosition(step)[2] << endl;
     241    *out << getPositionAtStep(step)[0] << "\t";
     242    *out << getPositionAtStep(step)[1] << "\t";
     243    *out << getPositionAtStep(step)[2] << endl;
    244244    return true;
    245245  } else

  • src/atom_atominfo.cpp

    r6625c3 r056e70  
    8080const double& AtomInfo::operator[](size_t i) const
    8181{
     82  ASSERT(AtomicPosition.size() > 0,
     83      "AtomInfo::operator[]() - Access out of range.");
    8284  return AtomicPosition[0][i];
    8385}
    … …  
    8587const double& AtomInfo::at(size_t i) const
    8688{
     89  ASSERT(AtomicPosition.size() > 0,
     90      "AtomInfo::at() - Access out of range.");
    8791  return AtomicPosition[0].at(i);
    8892}
    8993
     94const double& AtomInfo::atStep(size_t i, int _step) const
     95{
     96  ASSERT(AtomicPosition.size() > _step,
     97      "AtomInfo::atStep() - Access out of range.");
     98  return AtomicPosition[_step].at(i);
     99}
     100
    90101void AtomInfo::set(size_t i, const double value)
    91102{
     103  ASSERT(AtomicPosition.size() > 0,
     104      "AtomInfo::set() - Access out of range.");
    92105  AtomicPosition[0].at(i) = value;
    93106}
    … …  
    95108const Vector& AtomInfo::getPosition() const
    96109{
     110  ASSERT(AtomicPosition.size() > 0,
     111      "AtomInfo::getPosition() - Access out of range.");
    97112  return AtomicPosition[0];
    98113}
    99114
    100 const Vector& AtomInfo::getPosition(const int _step) const
     115const Vector& AtomInfo::getPositionAtStep(const int _step) const
    101116{
    102117  ASSERT(_step < AtomicPosition.size(),
    103       "AtomInfo::getPosition() - Access out of range.");
     118      "AtomInfo::getPositionAtStep() - Access out of range.");
    104119  return AtomicPosition[_step];
    105120}
    … …  
    114129}
    115130
    116 Vector& AtomInfo::getAtomicVelocity()
    117 {
     131//Vector& AtomInfo::getAtomicVelocity()
     132//{
     133//  return AtomicVelocity[0];
     134//}
     135
     136//Vector& AtomInfo::getAtomicVelocity(const int _step)
     137//{
     138//  ASSERT(_step < AtomicVelocity.size(),
     139//      "AtomInfo::getAtomicVelocity() - Access out of range.");
     140//  return AtomicVelocity[_step];
     141//}
     142
     143const Vector& AtomInfo::getAtomicVelocity() const
     144{
     145  ASSERT(AtomicVelocity.size() > 0,
     146      "AtomInfo::getAtomicVelocity() - Access out of range.");
    118147  return AtomicVelocity[0];
    119148}
    120149
    121 Vector& AtomInfo::getAtomicVelocity(const int _step)
     150const Vector& AtomInfo::getAtomicVelocityAtStep(const int _step) const
    122151{
    123152  ASSERT(_step < AtomicVelocity.size(),
    … …  
    126155}
    127156
    128 const Vector& AtomInfo::getAtomicVelocity() const
    129 {
    130   return AtomicVelocity[0];
    131 }
    132 
    133 const Vector& AtomInfo::getAtomicVelocity(const int _step) const
    134 {
    135   ASSERT(_step < AtomicVelocity.size(),
    136       "AtomInfo::getAtomicVelocity() - Access out of range.");
    137   return AtomicVelocity[_step];
    138 }
    139 
    140157void AtomInfo::setAtomicVelocity(const Vector &_newvelocity)
    141158{
     159  ASSERT(0 < AtomicVelocity.size(),
     160      "AtomInfo::setAtomicVelocity() - Access out of range.");
    142161  AtomicVelocity[0] = _newvelocity;
    143162}
    144163
    145 void AtomInfo::setAtomicVelocity(const int _step, const Vector &_newvelocity)
     164void AtomInfo::setAtomicVelocityAtStep(const int _step, const Vector &_newvelocity)
    146165{
    147166  ASSERT(_step <= AtomicVelocity.size(),
    148       "AtomInfo::setAtomicVelocity() - Access out of range.");
     167      "AtomInfo::setAtomicVelocityAtStep() - Access out of range.");
    149168  if(_step < (int)AtomicVelocity.size()) {
    150169    AtomicVelocity[_step] = _newvelocity;
    … …  
    156175const Vector& AtomInfo::getAtomicForce() const
    157176{
     177  ASSERT(0 < AtomicForce.size(),
     178      "AtomInfo::getAtomicForce() - Access out of range.");
    158179  return AtomicForce[0];
    159180}
    160181
    161 const Vector& AtomInfo::getAtomicForce(const int _step) const
     182const Vector& AtomInfo::getAtomicForceAtStep(const int _step) const
    162183{
    163184  ASSERT(_step < AtomicForce.size(),
    … …  
    168189void AtomInfo::setAtomicForce(const Vector &_newforce)
    169190{
     191  ASSERT(0 < AtomicForce.size(),
     192      "AtomInfo::setAtomicForce() - Access out of range.");
    170193  AtomicForce[0] = _newforce;
    171194}
    172195
    173 void AtomInfo::setAtomicForce(const int _step, const Vector &_newforce)
    174 {
    175   ASSERT(_step <= AtomicForce.size(),
     196void AtomInfo::setAtomicForceAtStep(const int _step, const Vector &_newforce)
     197{
     198  const int size = AtomicForce.size();
     199  ASSERT(_step <= size,
    176200      "AtomInfo::setAtomicForce() - Access out of range.");
    177   if(_step < (int)AtomicForce.size()) {
     201  if(_step < size) {
    178202    AtomicForce[_step] = _newforce;
    179   } else if (_step == (int)AtomicForce.size()) {
     203  } else if (_step == size) {
    180204    AtomicForce.push_back(_newforce);
    181205  }
    … …  
    194218void AtomInfo::setPosition(const Vector& _vector)
    195219{
     220  ASSERT(0 < AtomicPosition.size(),
     221      "AtomInfo::setPosition() - Access out of range.");
    196222  AtomicPosition[0] = _vector;
    197223  //cout << "AtomInfo::setPosition: " << getType()->symbol << " at " << getPosition() << endl;
    198224}
    199225
    200 void AtomInfo::setPosition(int _step, const Vector& _vector)
     226void AtomInfo::setPositionAtStep(int _step, const Vector& _vector)
    201227{
    202228  ASSERT(_step <= AtomicPosition.size(),
    … …  
    212238const VectorInterface& AtomInfo::operator+=(const Vector& b)
    213239{
     240  ASSERT(0 < AtomicPosition.size(),
     241      "AtomInfo::operator+=() - Access out of range.");
    214242  AtomicPosition[0] += b;
    215243  return *this;
    … …  
    218246const VectorInterface& AtomInfo::operator-=(const Vector& b)
    219247{
     248  ASSERT(0 < AtomicPosition.size(),
     249      "AtomInfo::operator-=() - Access out of range.");
    220250  AtomicPosition[0] -= b;
    221251  return *this;
    … …  
    224254Vector const AtomInfo::operator+(const Vector& b) const
    225255{
     256  ASSERT(0 < AtomicPosition.size(),
     257      "AtomInfo::operator+() - Access out of range.");
    226258  Vector a(AtomicPosition[0]);
    227259  a += b;
    … …  
    231263Vector const AtomInfo::operator-(const Vector& b) const
    232264{
     265  ASSERT(0 < AtomicPosition.size(),
     266      "AtomInfo::operator-() - Access out of range.");
    233267  Vector a(AtomicPosition[0]);
    234268  a -= b;
    … …  
    238272double AtomInfo::distance(const Vector &point) const
    239273{
     274  ASSERT(0 < AtomicPosition.size(),
     275      "AtomInfo::distance() - Access out of range.");
    240276  return AtomicPosition[0].distance(point);
    241277}
    … …  
    243279double AtomInfo::DistanceSquared(const Vector &y) const
    244280{
     281  ASSERT(0 < AtomicPosition.size(),
     282      "AtomInfo::DistanceSquared() - Access out of range.");
    245283  return AtomicPosition[0].DistanceSquared(y);
    246284}
    … …  
    248286double AtomInfo::distance(const VectorInterface &_atom) const
    249287{
     288  ASSERT(0 < AtomicPosition.size(),
     289      "AtomInfo::distance() - Access out of range.");
    250290  return _atom.distance(AtomicPosition[0]);
    251291}
    … …  
    253293double AtomInfo::DistanceSquared(const VectorInterface &_atom) const
    254294{
     295  ASSERT(0 < AtomicPosition.size(),
     296      "AtomInfo::DistanceSquared() - Access out of range.");
    255297  return _atom.DistanceSquared(AtomicPosition[0]);
    256298}
    … …  
    258300VectorInterface &AtomInfo::operator=(const Vector& _vector)
    259301{
     302  ASSERT(0 < AtomicPosition.size(),
     303      "AtomInfo::operator=() - Access out of range.");
    260304  AtomicPosition[0] = _vector;
    261305  return *this;
    … …  
    264308void AtomInfo::ScaleAll(const double *factor)
    265309{
     310  ASSERT(0 < AtomicPosition.size(),
     311      "AtomInfo::ScaleAll() - Access out of range.");
    266312  AtomicPosition[0].ScaleAll(factor);
    267313}
    … …  
    269315void AtomInfo::ScaleAll(const Vector &factor)
    270316{
     317  ASSERT(0 < AtomicPosition.size(),
     318      "AtomInfo::ScaleAll() - Access out of range.");
    271319  AtomicPosition[0].ScaleAll(factor);
    272320}
    … …  
    274322void AtomInfo::Scale(const double factor)
    275323{
     324  ASSERT(0 < AtomicPosition.size(),
     325      "AtomInfo::Scale() - Access out of range.");
    276326  AtomicPosition[0].Scale(factor);
    277327}
    … …  
    279329void AtomInfo::Zero()
    280330{
     331  ASSERT(0 < AtomicPosition.size(),
     332      "AtomInfo::Zero() - Access out of range.");
    281333  AtomicPosition[0].Zero();
    282334}
    … …  
    284336void AtomInfo::One(const double one)
    285337{
     338  ASSERT(0 < AtomicPosition.size(),
     339      "AtomInfo::One() - Access out of range.");
    286340  AtomicPosition[0].One(one);
    287341}
    … …  
    289343void AtomInfo::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
    290344{
     345  ASSERT(0 < AtomicPosition.size(),
     346      "AtomInfo::LinearCombinationOfVectors() - Access out of range.");
    291347  AtomicPosition[0].LinearCombinationOfVectors(x1,x2,x3,factors);
    292348}
    … …  
    295351 *  returns the kinetic energy of this atom at a given time step
    296352 */
    297 double AtomInfo::getKineticEnergy(unsigned int step) const{
    298   return getMass() * AtomicVelocity[step].NormSquared();
    299 }
    300 
    301 Vector AtomInfo::getMomentum(unsigned int step) const{
    302   return getMass()*AtomicVelocity[step];
     353double AtomInfo::getKineticEnergy(unsigned int _step) const{
     354  ASSERT(_step < AtomicPosition.size(),
     355      "AtomInfo::getKineticEnergy() - Access out of range.");
     356  return getMass() * AtomicVelocity[_step].NormSquared();
     357}
     358
     359Vector AtomInfo::getMomentum(unsigned int _step) const{
     360  ASSERT(_step < AtomicPosition.size(),
     361      "AtomInfo::getMomentum() - Access out of range.");
     362  return getMass()*AtomicVelocity[_step];
    303363}
    304364
    … …  
    354414void AtomInfo::VelocityVerletUpdate(int nr, int NextStep, config *configuration, ForceMatrix *Force, const size_t offset)
    355415{
    356   //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
    357   if ((int)AtomicPosition.size() <= NextStep)
    358     ResizeTrajectory(NextStep+10);
    359   for (int d=0; d<NDIM; d++) {
    360     AtomicForce.at(NextStep)[d] = -Force->Matrix[0][nr][d+offset]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
    361     AtomicPosition.at(NextStep)[d] = AtomicPosition.at(NextStep-1)[d];
    362     AtomicPosition.at(NextStep)[d] += configuration->Deltat*(AtomicVelocity.at(NextStep-1)[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
    363     AtomicPosition.at(NextStep)[d] += 0.5*configuration->Deltat*configuration->Deltat*(AtomicForce.at(NextStep)[d]/getMass());     // F = m * a and s =
    364   }
     416  // update force
     417  // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
     418  Vector tempVector;
     419  for (int d=0; d<NDIM; d++)
     420    tempVector[d] = -Force->Matrix[0][nr][d+offset]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
     421  setAtomicForceAtStep(NextStep, tempVector);
     422
     423  // update position
     424  tempVector = getPositionAtStep(NextStep-1);
     425  tempVector += configuration->Deltat*(getAtomicVelocityAtStep(NextStep-1));     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
     426  tempVector += 0.5*configuration->Deltat*configuration->Deltat*(getAtomicForceAtStep(NextStep))*(1./getMass());     // F = m * a and s =
     427  setPositionAtStep(NextStep, tempVector);
     428
    365429  // Update U
    366   for (int d=0; d<NDIM; d++) {
    367     AtomicVelocity.at(NextStep)[d] = AtomicVelocity.at(NextStep-1)[d];
    368     AtomicVelocity.at(NextStep)[d] += configuration->Deltat * (AtomicForce.at(NextStep)[d]+AtomicForce.at(NextStep-1)[d]/getMass()); // v = F/m * t
    369   }
    370   // Update R (and F)
    371 //      out << "Integrated position&velocity of step " << (NextStep) << ": (";
    372 //      for (int d=0;d<NDIM;d++)
    373 //        out << AtomicPosition.at(NextStep).x[d] << " ";          // next step
    374 //      out << ")\t(";
    375 //      for (int d=0;d<NDIM;d++)
    376 //        Log() << Verbose(0) << AtomicVelocity.at(NextStep).x[d] << " ";          // next step
    377 //      out << ")" << endl;
     430  tempVector = getAtomicVelocityAtStep(NextStep-1);
     431  tempVector += configuration->Deltat * (getAtomicForceAtStep(NextStep)+getAtomicForceAtStep(NextStep-1))*(1./getMass()); // v = F/m * t
     432  setAtomicVelocityAtStep(NextStep, tempVector);
     433
     434  // some info for debugging
     435  DoLog(2) && (Log() << Verbose(2)
     436      << "Integrated position&velocity of step " << (NextStep) << ": ("
     437      << getPositionAtStep(NextStep) << ")\t("
     438      << getAtomicVelocityAtStep(NextStep) << ")" << std::endl);
    378439};
    379440

  • src/atom_atominfo.hpp

    r6625c3 r056e70  
    6464   * @return constant reference to AtomicVelocity
    6565   */
    66   Vector& getAtomicVelocity();
    67   /** Getter for AtomicVelocity.
    68    *
    69    * @param _step time step to return
    70    * @return constant reference to AtomicVelocity
    71    */
    72   Vector& getAtomicVelocity(const int _step);
     66//  Vector& getAtomicVelocity();
     67  /** Getter for AtomicVelocity.
     68   *
     69   * @param _step time step to return
     70   * @return constant reference to AtomicVelocity
     71   */
     72//  Vector& getAtomicVelocity(const int _step);
    7373  /** Getter for AtomicVelocity.
    7474   *
    … …  
    8383   * @return constant reference to AtomicVelocity
    8484   */
    85   const Vector& getAtomicVelocity(const int _step) const;
     85  const Vector& getAtomicVelocityAtStep(const int _step) const;
    8686  /** Setter for AtomicVelocity.
    8787   *
    … …  
    9696   * @param _newvelocity new velocity to set
    9797   */
    98   void setAtomicVelocity(const int _step, const Vector &_newvelocity);
     98  void setAtomicVelocityAtStep(const int _step, const Vector &_newvelocity);
    9999
    100100  /** Getter for AtomicForce.
    … …  
    110110   * @return constant reference to AtomicForce
    111111   */
    112   const Vector& getAtomicForce(const int _step) const;
     112  const Vector& getAtomicForceAtStep(const int _step) const;
    113113  /** Setter for AtomicForce.
    114114   *
    … …  
    123123   * @param _newvelocity new force vector to set
    124124   */
    125   void setAtomicForce(const int _step, const Vector &_newforce);
     125  void setAtomicForceAtStep(const int _step, const Vector &_newforce);
    126126
    127127  /** Getter for FixedIon.
    … …  
    165165   * @return atomic position at time step _step
    166166   */
    167   const Vector& atStep(size_t i, int _step) const;
     167  const double& atStep(size_t i, int _step) const;
    168168  /** Setter for AtomicPosition.
    169169   *
    … …  
    193193   * @return atomic position at time step _step
    194194   */
    195   const Vector& getPosition(int _step) const;
     195  const Vector& getPositionAtStep(int _step) const;
    196196
    197197  // Assignment operator
    … …  
    208208   * @param _vector new position to set for time step _step
    209209   */
    210   void setPosition(int _step, const Vector& _vector);
     210  void setPositionAtStep(int _step, const Vector& _vector);
    211211  class VectorInterface &operator=(const Vector& _vector);
    212212

  • src/config.cpp

    r6625c3 r056e70  
    530530            if (!status)
    531531              break;
    532             neues->setFixedIon(_fixedion);
    533             neues->setPosition(position);
    534532
    535533            // check size of vectors
    … …  
    540538              neues->ResizeTrajectory(repetition+1);
    541539            }
     540            neues->setFixedIon(_fixedion);
     541            neues->setPositionAtStep(repetition,position);
    542542
    543543            // parse velocities if present
    … …  
    549549            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &tempVector[2], 1,optional))
    550550              tempVector[2] = 0.;
    551             neues->setAtomicVelocity(repetition, tempVector);
     551            neues->setAtomicVelocityAtStep(repetition, tempVector);
    552552
    553553            // parse forces if present
    … …  
    559559            if(!ParseForParameter(verbose,FileBuffer, keyword, 1, 10, 1, double_type, &tempVector[2], 1,optional))
    560560              tempVector[2] = 0.;
    561             neues->setAtomicForce(repetition, tempVector);
     561            neues->setAtomicForceAtStep(repetition, tempVector);
    562562
    563563  //            Log() << Verbose(0) << "Parsed position of step " << (repetition) << ": (";

  • src/molecule.cpp

    r6625c3 r056e70  
    618618      Walker->setType(1);
    619619    }
    620     if (Walker->getTrajectorySize() <= (unsigned int)MDSteps) {
    621       Walker->ResizeTrajectory(MDSteps+10);
    622     }
     620
    623621    Walker->setPosition(Vector(x));
    624     Walker->setPosition(MDSteps-1, Vector(x));
    625     Walker->setAtomicVelocity(MDSteps-1, zeroVec);
    626     Walker->setAtomicForce(MDSteps-1, zeroVec);
     622    Walker->setPositionAtStep(MDSteps-1, Vector(x));
     623    Walker->setAtomicVelocityAtStep(MDSteps-1, zeroVec);
     624    Walker->setAtomicForceAtStep(MDSteps-1, zeroVec);
    627625    AddAtom(Walker);  // add to molecule
    628626    delete(item);

  • src/molecule_dynamics.cpp

    r6625c3 r056e70  
    6060    // determine normalized trajectories direction vector (n1, n2)
    6161    Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
    62     trajectory1 = Sprinter->getPosition(Params.endstep) - Walker->getPosition(Params.startstep);
     62    trajectory1 = Sprinter->getPositionAtStep(Params.endstep) - Walker->getPositionAtStep(Params.startstep);
    6363    trajectory1.Normalize();
    6464    Norm1 = trajectory1.Norm();
    6565    Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
    66     trajectory2 = Sprinter->getPosition(Params.endstep) - (*iter)->getPosition(Params.startstep);
     66    trajectory2 = Sprinter->getPositionAtStep(Params.endstep) - (*iter)->getPositionAtStep(Params.startstep);
    6767    trajectory2.Normalize();
    6868    Norm2 = trajectory1.Norm();
    6969    // check whether either is zero()
    7070    if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
    71       tmp = Walker->getPosition(Params.startstep).distance((*iter)->getPosition(Params.startstep));
     71      tmp = Walker->getPositionAtStep(Params.startstep).distance((*iter)->getPositionAtStep(Params.startstep));
    7272    } else if (Norm1 < MYEPSILON) {
    7373      Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
    74       trajectory1 = Sprinter->getPosition(Params.endstep) - (*iter)->getPosition(Params.startstep);
     74      trajectory1 = Sprinter->getPositionAtStep(Params.endstep) - (*iter)->getPositionAtStep(Params.startstep);
    7575      trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    7676      trajectory1 -= trajectory2;   // project the part in norm direction away
    … …  
    7878    } else if (Norm2 < MYEPSILON) {
    7979      Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
    80       trajectory2 = Sprinter->getPosition(Params.endstep) - Walker->getPosition(Params.startstep);  // copy second offset
     80      trajectory2 = Sprinter->getPositionAtStep(Params.endstep) - Walker->getPositionAtStep(Params.startstep);  // copy second offset
    8181      trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
    8282      trajectory2 -= trajectory1;   // project the part in norm direction away
    … …  
    8787  //        Log() << Verbose(0) << " and ";
    8888  //        Log() << Verbose(0) << trajectory2;
    89       tmp = Walker->getPosition(Params.startstep).distance((*iter)->getPosition(Params.startstep));
     89      tmp = Walker->getPositionAtStep(Params.startstep).distance((*iter)->getPositionAtStep(Params.startstep));
    9090  //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
    9191    } else { // determine distance by finding minimum distance
    … …  
    106106        gsl_matrix_set(A, 1, i, trajectory2[i]);
    107107        gsl_matrix_set(A, 2, i, normal[i]);
    108         gsl_vector_set(x,i, (Walker->getPosition(Params.startstep)[i] - (*iter)->getPosition(Params.startstep)[i]));
     108        gsl_vector_set(x,i, (Walker->getPositionAtStep(Params.startstep)[i] - (*iter)->getPositionAtStep(Params.startstep)[i]));
    109109      }
    110110      // solve the linear system by Householder transformations
    … …  
    117117      trajectory2.Scale(gsl_vector_get(x,1));
    118118      normal.Scale(gsl_vector_get(x,2));
    119       TestVector = (*iter)->getPosition(Params.startstep) + trajectory2 + normal
    120                    - (Walker->getPosition(Params.startstep) + trajectory1);
     119      TestVector = (*iter)->getPositionAtStep(Params.startstep) + trajectory2 + normal
     120                   - (Walker->getPositionAtStep(Params.startstep) + trajectory1);
    121121      if (TestVector.Norm() < MYEPSILON) {
    122122  //          Log() << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;
    … …  
    183183    // first term: distance to target
    184184    Runner = Params.PermutationMap[(*iter)->nr];   // find target point
    185     tmp = ((*iter)->getPosition(Params.startstep).distance(Runner->getPosition(Params.endstep)));
     185    tmp = ((*iter)->getPositionAtStep(Params.startstep).distance(Runner->getPositionAtStep(Params.endstep)));
    186186    tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    187187    result += Params.PenaltyConstants[0] * tmp;
    … …  
    239239  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    240240    for (molecule::const_iterator runner = mol->begin(); runner != mol->end(); ++runner) {
    241       Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->getPosition(Params.startstep).distance((*runner)->getPosition(Params.endstep)), (*runner)) );
     241      Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->getPositionAtStep(Params.startstep).distance((*runner)->getPositionAtStep(Params.endstep)), (*runner)) );
    242242    }
    243243  }
    … …  
    489489    // set forces
    490490    for (int i=NDIM;i++;)
    491       Force->Matrix[0][_atom->nr][5+i] += 2.*constant*sqrt(_atom->getPosition(startstep).distance(Sprinter->getPosition(endstep)));
     491      Force->Matrix[0][_atom->nr][5+i] += 2.*constant*sqrt(_atom->getPositionAtStep(startstep).distance(Sprinter->getPositionAtStep(endstep)));
    492492  }
    493493  DoLog(1) && (Log() << Verbose(1) << "done." << endl);
    … …  
    540540      Sprinter = mol->AddCopyAtom((*iter));
    541541      // add to Trajectories
    542       Vector temp = (*iter)->getPosition(startstep) + (PermutationMap[(*iter)->nr]->getPosition(endstep) - (*iter)->getPosition(startstep))*((double)step/(double)MaxSteps);
     542      Vector temp = (*iter)->getPositionAtStep(startstep) + (PermutationMap[(*iter)->nr]->getPositionAtStep(endstep) - (*iter)->getPositionAtStep(startstep))*((double)step/(double)MaxSteps);
    543543      Sprinter->setPosition(temp);
    544       (*iter)->setAtomicVelocity(step, zeroVec);
    545       (*iter)->setAtomicForce(step, zeroVec);
     544      (*iter)->setAtomicVelocityAtStep(step, zeroVec);
     545      (*iter)->setAtomicForceAtStep(step, zeroVec);
    546546      //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
    547547    }

  • src/molecule_geometry.cpp

    r6625c3 r056e70  
    317317  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    318318    for (size_t j=0;j<(*iter)->getTrajectorySize();j++) {
    319       Vector temp = (*iter)->getPosition(j);
     319      Vector temp = (*iter)->getPositionAtStep(j);
    320320      temp.ScaleAll(*factor);
    321       (*iter)->setPosition(j,temp);
     321      (*iter)->setPositionAtStep(j,temp);
    322322    }
    323323  }
    … …  
    331331  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    332332    for (size_t j=0;j<(*iter)->getTrajectorySize();j++) {
    333       Vector temp = (*iter)->getPosition(j);
    334       temp += (*trans);
    335       (*iter)->setPosition(j, temp);
     333      (*iter)->setPositionAtStep(j, (*iter)->getPositionAtStep(j) + (*trans));
    336334    }
    337335  }
    … …  
    445443    for (int j=0;j<MDSteps;j++) {
    446444      Vector temp;
    447       temp[0] =  cos(alpha) * (*iter)->getPosition(j)[0] + sin(alpha) * (*iter)->getPosition(j)[2];
    448       temp[2] = -sin(alpha) * (*iter)->getPosition(j)[0] + cos(alpha) * (*iter)->getPosition(j)[2];
    449       (*iter)->setPosition(j,temp);
     445      temp[0] =  cos(alpha) * (*iter)->getPositionAtStep(j)[0] + sin(alpha) * (*iter)->getPositionAtStep(j)[2];
     446      temp[2] = -sin(alpha) * (*iter)->getPositionAtStep(j)[0] + cos(alpha) * (*iter)->getPositionAtStep(j)[2];
     447      (*iter)->setPositionAtStep(j,temp);
    450448    }
    451449  }
    … …  
    465463    for (int j=0;j<MDSteps;j++) {
    466464      Vector temp;
    467       temp[1] =  cos(alpha) * (*iter)->getPosition(j)[1] + sin(alpha) * (*iter)->getPosition(j)[2];
    468       temp[2] = -sin(alpha) * (*iter)->getPosition(j)[1] + cos(alpha) * (*iter)->getPosition(j)[2];
    469       (*iter)->setPosition(j,temp);
     465      temp[1] =  cos(alpha) * (*iter)->getPositionAtStep(j)[1] + sin(alpha) * (*iter)->getPositionAtStep(j)[2];
     466      temp[2] = -sin(alpha) * (*iter)->getPositionAtStep(j)[1] + cos(alpha) * (*iter)->getPositionAtStep(j)[2];
     467      (*iter)->setPositionAtStep(j,temp);
    470468    }
    471469  }
Note: See TracChangeset for help on using the changeset viewer.