Changes in / [c7fe90:5c8807]


Ignore:
Files:
18 added
30 edited

Legend:

Unmodified
Added
Removed

  • doc/userguide/userguide.xml

    rc7fe90 r5c8807  
    955955          domain.</para>
    956956        </section>
     957
     958        <section xml:id='atoms.mirror-atoms'>
     959          <title xml:id='atoms.mirror-atoms.title'>Mirroring atoms</title>
     960
     961          <para>Present (and selected) atoms can be mirrored with respect to
     962          a certain plane. You have to specify the normal vector of the plane
     963          and the offset with respect to the origin as follows</para>
     964
     965          <programlisting>
     966                ... --mirror-atoms "1,0,0" \
     967                    --plane-offset 10.1 \
     968                    --periodic 0
     969          </programlisting>
     970        </section>
    957971
    958972        <section xml:id='atoms.translate-to-origin'>
    … …  
    21332147            a fragment of order 1, e.g. a single hydrogen atom.</para>
    21342148          </note>
     2149        </section>
     2150
     2151        <section xml:id='potentials.fit-compound-potential'>
     2152          <title xml:id='potentials.fit-compound-potential.title'>Fitting
     2153          many empirical potentials simultaneously</title>
     2154
    21352155
    21362156          <para>Another way is using a file containing a specific set of
    … …  
    21382158
    21392159          <programlisting>
    2140                 ... --fit-potential \
     2160                ... --fit-compound-potential \
    21412161                    --fragment-charges 8 1 1 \
    21422162                    --potential-file water.potentials \
    … …  
    21602180          type of analysis.</para>
    21612181
    2162           <para>Note that you can combine the two ways, i.e. start with the
    2163           first but give an empty potential file. The resulting parameters are
    2164           stored in this way. Fit other potentials and give different file
    2165           names for each. Eventually, you have to combine the file in a text
    2166           editor at the moment.</para>
     2182          <para>Note that you can combine the two ways, i.e. start with a
     2183          fit-potential call but give an empty potential file. The resulting
     2184          parameters are stored in it. Fit other potentials and give different
     2185          file names for each in turn. Eventually, you have to combine the file
     2186          in a text editor at the moment. And perform a fit-compound-potential
     2187          with this file.</para>
     2188        </section>
     2189
     2190       
     2191        <section xml:id='potentials.parse-potential'>
     2192          <title xml:id='potentials.parse-potential.title'>Parsing an
     2193          empirical potentials file</title>
     2194
     2195          <para>Responsible for the compound potential is every potential
     2196          function whose signature matches with the designated fragment-charges
     2197          and who is currently known to an internal instance called the
     2198          PotentialRegistry.</para>
     2199
     2200          <para>More potentials can be registered (fit-potential will also
     2201          register the potential it fits) by parsing them from a file.</para>
     2202
     2203          <programlisting>
     2204                ... --parse-potentials water.potentials
     2205          </programlisting>
     2206
     2207          <note>Currently, only <productname>TREMOLO</productname> potential
     2208          files are understood and can be parsed.</note>
     2209        </section>
     2210
     2211        <section xml:id='potentials.save-potential'>
     2212          <title xml:id='potentials.save-potential.title'>Saving an
     2213          empirical potentials file</title>
     2214
     2215          <para>The opposite to parse-potentials is save-potentials that writes
     2216          every potential currently known to the PotentialRegistry to the given
     2217          file along with the currently fitted parameters</para>
     2218
     2219          <programlisting>
     2220                ... --save-potentials water.potentials
     2221          </programlisting>
     2222
     2223          <note>Again, only the <productname>TREMOLO</productname> potential
     2224          format is understood currently and is written.</note>
    21672225        </section>
    21682226

  • src/Actions/GlobalListOfActions.hpp

    rc7fe90 r5c8807  
    3232  (AtomAdd) \
    3333  (AtomChangeElement) \
     34  (AtomMirror) \
    3435  (AtomRemove) \
    3536  (AtomRotateAroundOriginByAngle) \
    … …  
    8586  (PotentialFitParticleCharges) \
    8687  (PotentialParseHomologies) \
     88  (PotentialParsePotentials) \
    8789  (PotentialSaveHomologies) \
     90  (PotentialSavePotentials) \
    8891  (ParserParseTremoloPotentials) \
    8992  (ParserSaveSelectedAtomsAsExtTypes) \
    … …  
    150153#define GLOBALLISTOFACTIONS_LEVMAR \
    151154    BOOST_PP_SEQ_PUSH_BACK( \
    152         GLOBALLISTOFACTIONS_initial, \
    153         PotentialFitPotential \
     155                BOOST_PP_SEQ_PUSH_BACK( \
     156                                GLOBALLISTOFACTIONS_initial, \
     157                                PotentialFitPotential \
     158                ), \
     159        PotentialFitCompoundPotential \
    154160    )
    155161#else

  • src/Actions/Makefile.am

    rc7fe90 r5c8807  
    147147  Actions/AtomAction/AddAction.cpp \
    148148  Actions/AtomAction/ChangeElementAction.cpp \
     149  Actions/AtomAction/MirrorAction.cpp \
    149150  Actions/AtomAction/RemoveAction.cpp \
    150151  Actions/AtomAction/RotateAroundOriginByAngleAction.cpp \
    … …  
    155156  Actions/AtomAction/AddAction.hpp \
    156157  Actions/AtomAction/ChangeElementAction.hpp \
     158  Actions/AtomAction/MirrorAction.hpp \
    157159  Actions/AtomAction/RemoveAction.hpp \
    158160  Actions/AtomAction/RotateAroundOriginByAngleAction.hpp \
    … …  
    163165  Actions/AtomAction/AddAction.def \
    164166  Actions/AtomAction/ChangeElementAction.def \
     167  Actions/AtomAction/MirrorAction.def \
    165168  Actions/AtomAction/RemoveAction.def \
    166169  Actions/AtomAction/RotateAroundOriginByAngleAction.def \
    … …  
    357360  Actions/PotentialAction/FitParticleChargesAction.cpp \
    358361  Actions/PotentialAction/ParseHomologiesAction.cpp \
    359   Actions/PotentialAction/SaveHomologiesAction.cpp
     362  Actions/PotentialAction/ParsePotentialsAction.cpp \
     363  Actions/PotentialAction/SaveHomologiesAction.cpp \
     364  Actions/PotentialAction/SavePotentialsAction.cpp
    360365POTENTIALACTIONHEADER = \
    361366  Actions/PotentialAction/FitParticleChargesAction.hpp \
    362367  Actions/PotentialAction/ParseHomologiesAction.hpp \
    363   Actions/PotentialAction/SaveHomologiesAction.hpp
     368  Actions/PotentialAction/ParsePotentialsAction.hpp \
     369  Actions/PotentialAction/SaveHomologiesAction.hpp \
     370  Actions/PotentialAction/SavePotentialsAction.hpp
    364371POTENTIALACTIONDEFS = \
    365372  Actions/PotentialAction/FitParticleChargesAction.def  \
    366373  Actions/PotentialAction/ParseHomologiesAction.def \
    367   Actions/PotentialAction/SaveHomologiesAction.def
     374  Actions/PotentialAction/ParsePotentialsAction.def \
     375  Actions/PotentialAction/SaveHomologiesAction.def \
     376  Actions/PotentialAction/SavePotentialsAction.def
    368377
    369378if CONDLEVMAR
    370379POTENTIALACTIONSOURCE += \
     380  Actions/PotentialAction/FitCompoundPotentialAction.cpp \
    371381  Actions/PotentialAction/FitPotentialAction.cpp
    372382POTENTIALACTIONHEADER += \
     383  Actions/PotentialAction/FitCompoundPotentialAction.hpp \
    373384  Actions/PotentialAction/FitPotentialAction.hpp
    374385POTENTIALACTIONDEFS += \
     386  Actions/PotentialAction/FitCompoundPotentialAction.def \
    375387  Actions/PotentialAction/FitPotentialAction.def
    376388endif

  • src/Actions/MoleculeAction/ForceAnnealingAction.cpp

    rc7fe90 r5c8807  
    7979        ++iter) {
    8080      atom * const Walker = iter->second;
    81       Walker->UpdateSteps();
    8281      Walker->setPositionAtStep(CurrentStep+1,
    8382          Walker->getPositionAtStep(CurrentStep));

  • src/Actions/MoleculeAction/VerletIntegrationAction.cpp

    rc7fe90 r5c8807  
    5656using namespace MoleCuilder;
    5757
    58 enum {
     58enum VectorIndexType {
    5959  PositionIndex=0,
    6060  VelocityIndex=1,
    6161  ForceIndex=2,
    6262  MAXINDEX
    63 } VectorIndexType;
     63};
    6464
    6565// and construct the stuff
    … …  
    134134
    135135  // remove current step for all modified atoms
    136   removeLastStep(getIdsFromAtomicInfo(state->UndoInfo));
     136  removeLastStep(getIdsFromAtomicInfo(state->UndoInfo), CurrentStep);
    137137
    138138  // and set back the old step (forces have been changed)
    … …  
    151151
    152152  // set stored new state
    153   addNewStep(state->UndoInfo);
     153  size_t CurrentStep = WorldTime::getInstance().getTime()+1;
     154  addNewStep(state->UndoInfo, CurrentStep);
    154155
    155156  // add a new time step
    156   size_t CurrentStep = WorldTime::getInstance().getTime();
    157   World::getInstance().setTime(CurrentStep+1);
     157  World::getInstance().setTime(CurrentStep);
    158158
    159159  // and set positions of the new step

  • src/Actions/PotentialAction/FitPotentialAction.cpp

    rc7fe90 r5c8807  
    5555#include "Fragmentation/Homology/HomologyGraph.hpp"
    5656#include "Fragmentation/Summation/SetValues/Fragment.hpp"
    57 #include "FunctionApproximation/Extractors.hpp"
    58 #include "FunctionApproximation/FunctionApproximation.hpp"
    59 #include "FunctionApproximation/FunctionModel.hpp"
    60 #include "FunctionApproximation/TrainingData.hpp"
    61 #include "FunctionApproximation/writeDistanceEnergyTable.hpp"
    62 #include "Potentials/CompoundPotential.hpp"
    63 #include "Potentials/Exceptions.hpp"
    64 #include "Potentials/PotentialDeserializer.hpp"
     57#include "Potentials/EmpiricalPotential.hpp"
    6558#include "Potentials/PotentialFactory.hpp"
    6659#include "Potentials/PotentialRegistry.hpp"
    6760#include "Potentials/PotentialSerializer.hpp"
     61#include "Potentials/PotentialTrainer.hpp"
    6862#include "Potentials/SerializablePotential.hpp"
     63#include "World.hpp"
    6964
    7065using namespace MoleCuilder;
    … …  
    7570/** =========== define the function ====================== */
    7671
    77 HomologyGraph getFirstGraphwithSpecifiedElements(
    78     const HomologyContainer &homologies,
    79     const SerializablePotential::ParticleTypes_t &types)
    80 {
    81   ASSERT( !types.empty(),
    82       "getFirstGraphwithSpecifiedElements() - charges is empty?");
    83   // create charges
    84   Fragment::charges_t charges;
    85   charges.resize(types.size());
    86   std::transform(types.begin(), types.end(),
    87       charges.begin(), boost::lambda::_1);
    88   // convert into count map
    89   Extractors::elementcounts_t counts_per_charge =
    90       Extractors::_detail::getElementCounts(charges);
    91   ASSERT( !counts_per_charge.empty(),
    92       "getFirstGraphwithSpecifiedElements() - charge counts are empty?");
    93   LOG(2, "DEBUG: counts_per_charge is " << counts_per_charge << ".");
    94   // we want to check each (unique) key only once
    95   for (HomologyContainer::const_key_iterator iter = homologies.key_begin();
    96       iter != homologies.key_end(); iter = homologies.getNextKey(iter)) {
    97     // check if every element has the right number of counts
    98     Extractors::elementcounts_t::const_iterator countiter = counts_per_charge.begin();
    99     for (; countiter != counts_per_charge.end(); ++countiter)
    100       if (!(*iter).hasTimesAtomicNumber(
    101           static_cast<size_t>(countiter->first),
    102           static_cast<size_t>(countiter->second))
    103           )
    104         break;
    105     if( countiter == counts_per_charge.end())
    106       return *iter;
    107   }
    108   return HomologyGraph();
    109 }
    110 
    111 SerializablePotential::ParticleTypes_t getNumbersFromElements(
    112     const std::vector<const element *> &fragment)
    113 {
    114   SerializablePotential::ParticleTypes_t fragmentnumbers;
    115   std::transform(fragment.begin(), fragment.end(), std::back_inserter(fragmentnumbers),
    116       boost::bind(&element::getAtomicNumber, _1));
    117   return fragmentnumbers;
    118 }
    119 
    120 
    12172ActionState::ptr PotentialFitPotentialAction::performCall() {
    12273  // fragment specifies the homology fragment to use
    12374  SerializablePotential::ParticleTypes_t fragmentnumbers =
    124       getNumbersFromElements(params.fragment.get());
     75      PotentialTrainer::getNumbersFromElements(params.fragment.get());
    12576
    12677  // either charges and a potential is specified or a file
    127   if (boost::filesystem::exists(params.potential_file.get())) {
    128     std::ifstream returnstream(params.potential_file.get().string().c_str());
    129     if (returnstream.good()) {
    130       try {
    131         PotentialDeserializer deserialize(returnstream);
    132         deserialize();
    133       } catch (SerializablePotentialMissingValueException &e) {
    134         if (const std::string *key = boost::get_error_info<SerializablePotentialKey>(e))
    135           STATUS("Missing value when parsing information for potential "+*key+".");
    136         else
    137           STATUS("Missing value parsing information for potential with unknown key.");
    138         return Action::failure;
    139       } catch (SerializablePotentialIllegalKeyException &e) {
    140         if (const std::string *key = boost::get_error_info<SerializablePotentialKey>(e))
    141           STATUS("Illegal key parsing information for potential "+*key+".");
    142         else
    143           STATUS("Illegal key parsing information for potential with unknown key.");
    144         return Action::failure;
    145       }
    146     } else {
    147       STATUS("Failed to parse from "+params.potential_file.get().string()+".");
    148       return Action::failure;
     78  if (params.charges.get().empty()) {
     79    STATUS("No charges given!");
     80    return Action::failure;
     81  } else {
     82    // charges specify the potential type
     83    SerializablePotential::ParticleTypes_t chargenumbers =
     84        PotentialTrainer::getNumbersFromElements(params.charges.get());
     85
     86    LOG(0, "STATUS: I'm training now a " << params.potentialtype.get()
     87        << " potential on charges " << chargenumbers << " on data from World's homologies.");
     88
     89    // register desired potential and an additional constant one
     90    {
     91      EmpiricalPotential *potential =
     92          PotentialFactory::getInstance().createInstance(
     93              params.potentialtype.get(),
     94              chargenumbers);
     95      // check whether such a potential already exists
     96      const std::string potential_name = potential->getName();
     97      if (PotentialRegistry::getInstance().isPresentByName(potential_name)) {
     98        delete potential;
     99        potential = PotentialRegistry::getInstance().getByName(potential_name);
     100      } else
     101        PotentialRegistry::getInstance().registerInstance(potential);
    149102    }
    150     returnstream.close();
    151 
    152     LOG(0, "STATUS: I'm training now a set of potentials parsed from "
    153         << params.potential_file.get().string() << " on a fragment "
    154         << fragmentnumbers << " on data from World's homologies.");
    155 
    156   } else {
    157     if (params.charges.get().empty()) {
    158       STATUS("Neither charges nor potential file given!");
    159       return Action::failure;
    160     } else {
    161       // charges specify the potential type
    162       SerializablePotential::ParticleTypes_t chargenumbers =
    163           getNumbersFromElements(params.charges.get());
    164 
    165       LOG(0, "STATUS: I'm training now a " << params.potentialtype.get()
    166           << " potential on charges " << chargenumbers << " on data from World's homologies.");
    167 
    168       // register desired potential and an additional constant one
    169       {
    170         EmpiricalPotential *potential =
    171             PotentialFactory::getInstance().createInstance(
    172                 params.potentialtype.get(),
    173                 chargenumbers);
    174         // check whether such a potential already exists
    175         const std::string potential_name = potential->getName();
    176         if (PotentialRegistry::getInstance().isPresentByName(potential_name)) {
    177           delete potential;
    178           potential = PotentialRegistry::getInstance().getByName(potential_name);
    179         } else
    180           PotentialRegistry::getInstance().registerInstance(potential);
    181       }
    182       {
    183         EmpiricalPotential *constant =
    184             PotentialFactory::getInstance().createInstance(
    185                 std::string("constant"),
    186                 SerializablePotential::ParticleTypes_t());
    187         // check whether such a potential already exists
    188         const std::string constant_name = constant->getName();
    189         if (PotentialRegistry::getInstance().isPresentByName(constant_name)) {
    190           delete constant;
    191           constant = PotentialRegistry::getInstance().getByName(constant_name);
    192         } else
    193           PotentialRegistry::getInstance().registerInstance(constant);
    194       }
     103    {
     104      EmpiricalPotential *constant =
     105          PotentialFactory::getInstance().createInstance(
     106              std::string("constant"),
     107              SerializablePotential::ParticleTypes_t());
     108      // check whether such a potential already exists
     109      const std::string constant_name = constant->getName();
     110      if (PotentialRegistry::getInstance().isPresentByName(constant_name)) {
     111        delete constant;
     112        constant = PotentialRegistry::getInstance().getByName(constant_name);
     113      } else
     114        PotentialRegistry::getInstance().registerInstance(constant);
    195115    }
    196116  }
    197117
    198118  // parse homologies into container
    199   HomologyContainer &homologies = World::getInstance().getHomologies();
     119  const HomologyContainer &homologies = World::getInstance().getHomologies();
    200120
    201121  // first we try to look into the HomologyContainer
    … …  
    211131
    212132  // then we ought to pick the right HomologyGraph ...
    213   const HomologyGraph graph = getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers);
     133  const HomologyGraph graph =
     134      PotentialTrainer::getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers);
    214135  if (graph != HomologyGraph()) {
    215136    LOG(1, "First representative graph containing fragment "
    … …  
    220141  }
    221142
    222   // fit potential
    223   FunctionModel *model = new CompoundPotential(graph);
    224   ASSERT( model != NULL,
    225       "PotentialFitPotentialAction::performCall() - model is NULL.");
    226 
    227   /******************** TRAINING ********************/
    228   // fit potential
    229   FunctionModel::parameters_t bestparams(model->getParameterDimension(), 0.);
    230   {
    231     // Afterwards we go through all of this type and gather the distance and the energy value
    232     TrainingData data(model->getSpecificFilter());
    233     data(homologies.getHomologousGraphs(graph));
    234 
    235     // print distances and energies if desired for debugging
    236     if (!data.getTrainingInputs().empty()) {
    237       // print which distance is which
    238       size_t counter=1;
    239       if (DoLog(3)) {
    240         const FunctionModel::arguments_t &inputs = data.getAllArguments()[0];
    241         for (FunctionModel::arguments_t::const_iterator iter = inputs.begin();
    242             iter != inputs.end(); ++iter) {
    243           const argument_t &arg = *iter;
    244           LOG(3, "DEBUG: distance " << counter++ << " is between (#"
    245               << arg.indices.first << "c" << arg.types.first << ","
    246               << arg.indices.second << "c" << arg.types.second << ").");
    247         }
    248       }
    249 
    250       // print table
    251       if (params.training_file.get().string().empty()) {
    252         LOG(3, "DEBUG: I gathered the following training data:\n" <<
    253             _detail::writeDistanceEnergyTable(data.getDistanceEnergyTable()));
    254       } else {
    255         std::ofstream trainingstream(params.training_file.get().string().c_str());
    256         if (trainingstream.good()) {
    257           LOG(3, "DEBUG: Writing training data to file " <<
    258               params.training_file.get().string() << ".");
    259           trainingstream << _detail::writeDistanceEnergyTable(data.getDistanceEnergyTable());
    260         }
    261         trainingstream.close();
    262       }
    263     }
    264 
    265     if ((params.threshold.get() < 1) && (params.best_of_howmany.isSet()))
    266       ELOG(2, "threshold parameter always overrules max_runs, both are specified.");
    267     // now perform the function approximation by optimizing the model function
    268     FunctionApproximation approximator(data, *model);
    269     if (model->isBoxConstraint() && approximator.checkParameterDerivatives()) {
    270       double l2error = std::numeric_limits<double>::max();
    271       // seed with current time
    272       srand((unsigned)time(0));
    273       unsigned int runs=0;
    274       // threshold overrules max_runs
    275       const double threshold = params.threshold.get();
    276       const unsigned int max_runs = (threshold >= 1.) ?
    277           (params.best_of_howmany.isSet() ? params.best_of_howmany.get() : 1) : 0;
    278       LOG(1, "INFO: Maximum runs is " << max_runs << " and threshold set to " << threshold << ".");
    279       do {
    280         // generate new random initial parameter values
    281         model->setParametersToRandomInitialValues(data);
    282         LOG(1, "INFO: Initial parameters of run " << runs << " are "
    283             << model->getParameters() << ".");
    284         approximator(FunctionApproximation::ParameterDerivative);
    285         LOG(1, "INFO: Final parameters of run " << runs << " are "
    286             << model->getParameters() << ".");
    287         const double new_l2error = data.getL2Error(*model);
    288         if (new_l2error < l2error) {
    289           // store currently best parameters
    290           l2error = new_l2error;
    291           bestparams = model->getParameters();
    292           LOG(1, "STATUS: New fit from run " << runs
    293               << " has better error of " << l2error << ".");
    294         }
    295       } while (( ++runs < max_runs) || (l2error > threshold));
    296       // reset parameters from best fit
    297       model->setParameters(bestparams);
    298       LOG(1, "INFO: Best parameters with L2 error of "
    299           << l2error << " are " << model->getParameters() << ".");
    300     } else {
    301       STATUS("No required parameter derivatives for a box constraint minimization known.");
    302       return Action::failure;
    303     }
    304 
    305     // create a map of each fragment with error.
    306     HomologyContainer::range_t fragmentrange = homologies.getHomologousGraphs(graph);
    307     TrainingData::L2ErrorConfigurationIndexMap_t WorseFragmentMap =
    308         data.getWorstFragmentMap(*model, fragmentrange);
    309     LOG(0, "RESULT: WorstFragmentMap " << WorseFragmentMap << ".");
    310 
    311     // print fitted potentials
    312     std::stringstream potentials;
    313     PotentialSerializer serialize(potentials);
    314     serialize();
    315     LOG(1, "STATUS: Resulting parameters are " << std::endl << potentials.str());
    316     std::ofstream returnstream(params.potential_file.get().string().c_str());
    317     if (returnstream.good()) {
    318       returnstream << potentials.str();
    319     }
     143  // training
     144  PotentialTrainer trainer;
     145  const bool status = trainer(
     146      homologies,
     147      graph,
     148      params.training_file.get(),
     149      params.threshold.get(),
     150      params.best_of_howmany.get());
     151  if (!status) {
     152    STATUS("No required parameter derivatives for a box constraint minimization known.");
     153    return Action::failure;
    320154  }
    321   delete model;
    322155
    323156  return Action::success;

  • src/Actions/PotentialAction/FitPotentialAction.def

    rc7fe90 r5c8807  
    1818#include "Parameters/Validators/Specific/ElementValidator.hpp"
    1919#include "Parameters/Validators/Specific/EmptyStringValidator.hpp"
    20 #include "Parameters/Validators/Specific/FileSuffixValidator.hpp"
    21 #include "Parameters/Validators/Specific/FilePresentValidator.hpp"
    2220#include "Parameters/Validators/Specific/PotentialTypeValidator.hpp"
    2321
    … …  
    2523// ValueStorage by the token "Z" -> first column: int, Z, "Z"
    2624// "undefine" if no parameters are required, use (NOPARAM_DEFAULT) for each (undefined) default value
    27 #define paramtypes (boost::filesystem::path)(std::string)(boost::filesystem::path)(std::vector<const element *>)(std::vector<const element *>)(unsigned int)(double)
    28 #define paramtokens ("training-file")("potential-type")("potential-file")("potential-charges")("fragment-charges")("take-best-of")("set-threshold")
    29 #define paramdescriptions ("optional file to write training data to")("potential type to fit")("optional potential file specifying multiple potentials to fit")("charges specifying the potential")("charges specifying the fragment")("take the best among this many approximations")("Require L2 error to be smaller than threshold, overrides number of attempts")
    30 #define paramdefaults (PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(PARAM_DEFAULT(3))(PARAM_DEFAULT(1.))
    31 #define paramreferences (training_file)(potentialtype)(potential_file)(charges)(fragment)(best_of_howmany)(threshold)
     25#define paramtypes (boost::filesystem::path)(std::string)(std::vector<const element *>)(std::vector<const element *>)(unsigned int)(double)
     26#define paramtokens ("training-file")("potential-type")("potential-charges")("fragment-charges")("take-best-of")("set-threshold")
     27#define paramdescriptions ("optional file to write training data to")("potential type to fit")("charges specifying the potential")("charges specifying the fragment")("take the best among this many approximations")("Require L2 error to be smaller than threshold, overrides number of attempts")
     28#define paramdefaults (PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(PARAM_DEFAULT(3))(PARAM_DEFAULT(1.))
     29#define paramreferences (training_file)(potentialtype)(charges)(fragment)(best_of_howmany)(threshold)
    3230#define paramvalids \
    3331(DummyValidator<boost::filesystem::path>()) \
    34 (EmptyStringValidator() || PotentialTypeValidator()) \
    35 (!FilePresentValidator() || FileSuffixValidator("potentials")) \
     32(PotentialTypeValidator()) \
    3633(STLVectorValidator< std::vector<const element *> >(0,99, ElementValidator())) \
    3734(STLVectorValidator< std::vector<const element *> >(1,99, ElementValidator())) \

  • src/Actions/SelectionAction/Atoms/AtomByIdAction.cpp

    rc7fe90 r5c8807  
    9393      break;
    9494    default:
    95       ASSERT(0, "SelectionAtomByIdAction::performCall() - this must not happen.");
     95      STATUS("No atoms have been selected.");
    9696      return Action::failure;
    9797      break;

  • src/Actions/SelectionAction/Atoms/NotAtomByIdAction.cpp

    rc7fe90 r5c8807  
    9393      break;
    9494    default:
    95       ASSERT(0, "SelectionAtomByIdAction::performCall() - this must not happen.");
     95      STATUS("No atoms have been selected.");
    9696      return Action::failure;
    9797      break;

  • src/Actions/UndoRedoHelpers.cpp

    rc7fe90 r5c8807  
    181181}
    182182
    183 void MoleCuilder::removeLastStep(const std::vector<atomId_t> &_atoms)
     183void MoleCuilder::removeLastStep(const std::vector<atomId_t> &_atoms, const unsigned int _step)
    184184{
    185185  for (size_t i=0; i<_atoms.size(); ++i) {
    186186    atom * const _atom = World::getInstance().getAtom(AtomById(_atoms[i]));
    187     _atom->removeSteps();
    188   }
    189 }
    190 
    191 void MoleCuilder::addNewStep(const std::vector<AtomicInfo> &_movedatoms)
     187    _atom->removeStep(_step);
     188  }
     189}
     190
     191void MoleCuilder::addNewStep(const std::vector<AtomicInfo> &_movedatoms, const unsigned int _step)
    192192{
    193193  for(size_t i=0; i< _movedatoms.size(); ++i) {
    194194    atom * const _atom = World::getInstance().getAtom(AtomById(_movedatoms[i].getId()));
    195     _atom->UpdateSteps();
    196   }
    197 }
    198 
    199 void MoleCuilder::addNewStep(const std::vector<atomId_t> &_ids)
     195    _atom->UpdateStep(_step);
     196  }
     197}
     198
     199void MoleCuilder::addNewStep(const std::vector<atomId_t> &_ids, const unsigned int _step)
    200200{
    201201  for(size_t i=0; i< _ids.size(); ++i) {
    202202    atom * const _atom = World::getInstance().getAtom(AtomById(_ids[i]));
    203     _atom->UpdateSteps();
     203    _atom->UpdateStep(_step);
    204204  }
    205205}

  • src/Actions/UndoRedoHelpers.hpp

    rc7fe90 r5c8807  
    107107   *
    108108   * @param movedatoms atoms whose last step in time to remove
     109   * @param _step which trajectory to remove
    109110   */
    110   void removeLastStep(const std::vector<atomId_t> &movedatoms);
     111  void removeLastStep(const std::vector<atomId_t> &movedatoms, const unsigned int _step);
    111112
    112113  /** Adds another time step to all \a movedatoms.
    … …  
    115116   *
    116117   * @param _ids atoms whose last step in time to remove
     118   * @param _step which trajectory to insert/assign
    117119   */
    118   void addNewStep(const std::vector<atomId_t> &_ids);
     120  void addNewStep(const std::vector<atomId_t> &_ids, const unsigned int _step);
    119121
    120122  /** Adds another time step to all \a movedatoms.
    … …  
    124126   *
    125127   * @param _movedatoms atoms whose last step in time to remove
     128   * @param _step which trajectory to insert/assign
    126129   */
    127   void addNewStep(const std::vector<AtomicInfo> &_movedatoms);
     130  void addNewStep(const std::vector<AtomicInfo> &_movedatoms, const unsigned int _step);
    128131
    129132  /** Helper function to extract id information from vector of AtomicInfo.

  • src/Atom/TesselPoint.cpp

    rc7fe90 r5c8807  
    5252{};
    5353
    54 void TesselPoint::UpdateSteps()
     54void TesselPoint::UpdateStep(const unsigned int _step)
    5555{
    56   ASSERT(0, "TesselPoint::UpdateSteps() should never be called, TesselPoints don't have trajectories.");
    57   AtomInfo::AppendTrajectoryStep();
     56  ASSERT(0, "TesselPoint::UpdateStep() should never be called, TesselPoints don't have trajectories.");
     57  AtomInfo::AppendTrajectoryStep(_step);
    5858};
    5959
    60 void TesselPoint::removeSteps()
     60void TesselPoint::removeStep(const unsigned int _step)
    6161{
    62   ASSERT(0, "TesselPoint::removeSteps() should never be called, TesselPoints don't have trajectories.");
    63   AtomInfo::removeTrajectoryStep();
     62  ASSERT(0, "TesselPoint::removeStep() should never be called, TesselPoints don't have trajectories.");
     63  AtomInfo::removeTrajectoryStep(_step);
    6464};
    6565

  • src/Atom/TesselPoint.hpp

    rc7fe90 r5c8807  
    3737   *
    3838   */
    39   virtual void UpdateSteps();
     39  virtual void UpdateStep(const unsigned int _step);
    4040
    4141  /** Pops the last step in all trajectory vectors.
    … …  
    4545   * the real functions, \sa removeTrajectoryStep(), by all necessary subclasses.
    4646   */
    47   virtual void removeSteps();
     47  virtual void removeStep(const unsigned int _step);
    4848
    4949  /** Getter for this.

  • src/Atom/atom.cpp

    rc7fe90 r5c8807  
    6969    mol(0)
    7070{
    71   AtomicPosition = pointer->AtomicPosition; // copy trajectory of coordination
    72   AtomicVelocity = pointer->AtomicVelocity; // copy trajectory of velocity
    73   AtomicForce = pointer->AtomicForce;
    7471  // sign on to global atom change tracker
    7572  AtomObserver::getInstance().AtomInserted(this);
    … …  
    9390
    9491
    95 void atom::UpdateSteps()
    96 {
    97   LOG(4,"atom::UpdateSteps() called.");
     92void atom::UpdateStep(const unsigned int _step)
     93{
     94  LOG(4,"atom::UpdateStep() called.");
    9895  // append to position, velocity and force vector
    99   AtomInfo::AppendTrajectoryStep();
     96  AtomInfo::AppendTrajectoryStep(WorldTime::getTime()+1);
    10097  // append to ListOfBonds vector
    101   BondedParticleInfo::AppendTrajectoryStep();
    102 }
    103 
    104 void atom::removeSteps()
    105 {
    106   LOG(4,"atom::removeSteps() called.");
     98  BondedParticleInfo::AppendTrajectoryStep(WorldTime::getTime()+1);
     99}
     100
     101void atom::removeStep(const unsigned int _step)
     102{
     103  LOG(4,"atom::removeStep() called.");
    107104  // append to position, velocity and force vector
    108   AtomInfo::removeTrajectoryStep();
     105  AtomInfo::removeTrajectoryStep(_step);
    109106  // append to ListOfBonds vector
    110   BondedParticleInfo::removeTrajectoryStep();
     107  BondedParticleInfo::removeTrajectoryStep(_step);
    111108}
    112109
    … …  
    313310atom* NewAtom(atomId_t _id){
    314311  atom * res = new atom();
    315   // extent trajectory to current time step
    316   const size_t CurrentTime = WorldTime::getTime();
    317   for (size_t step = res->getTrajectorySize(); step <= CurrentTime; ++step)
    318     res->UpdateSteps();
    319312  res->setId(_id);
    320313  return res;

  • src/Atom/atom.hpp

    rc7fe90 r5c8807  
    6565   * real functions, \sa AppendTrajectoryStep(), by all necessary subclasses.
    6666   */
    67   virtual void UpdateSteps();
     67  virtual void UpdateStep(const unsigned int _step);
    6868
    6969  /** Pops the last step in all trajectory vectors.
    … …  
    7373   * the real functions, \sa removeTrajectoryStep(), by all necessary subclasses.
    7474   */
    75   virtual void removeSteps();
     75  virtual void removeStep(const unsigned int _step);
    7676
    7777  /** Output of a single atom with given numbering.

  • src/Atom/atom_atominfo.cpp

    rc7fe90 r5c8807  
    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.
    56 *
    67 *
    … …  
    5556  charge(0.)
    5657{
    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 };
     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}
    6562
    6663/** Copy constructor of class AtomInfo.
    … …  
    6865AtomInfo::AtomInfo(const AtomInfo &_atom) :
    6966    AtomicPosition(_atom.AtomicPosition),
     67    AtomicVelocity(_atom.AtomicVelocity),
     68    AtomicForce(_atom.AtomicForce),
    7069    AtomicElement(_atom.AtomicElement),
    7170    FixedIon(_atom.FixedIon),
    7271    charge(_atom.charge)
    7372{
    74   AtomicVelocity.reserve(1);
    75   AtomicVelocity.push_back(zeroVec);
    76   AtomicForce.reserve(1);
    77   AtomicForce.push_back(zeroVec);
    78 };
     73}
    7974
    8075AtomInfo::AtomInfo(const VectorInterface &_v) :
    … …  
    8378    charge(0.)
    8479{
    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);
     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) );
    9283};
    9384
    … …  
    9889};
    9990
    100 void AtomInfo::AppendTrajectoryStep()
    101 {
    102   NOTIFY(TrajectoryChanged);
    103   AtomicPosition.push_back(zeroVec);
    104   AtomicVelocity.push_back(zeroVec);
    105   AtomicForce.push_back(zeroVec);
     91void 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) );
    10698  LOG(5,"AtomInfo::AppendTrajectoryStep() called, size is ("
    10799      << AtomicPosition.size() << ","
    … …  
    110102}
    111103
    112 void AtomInfo::removeTrajectoryStep()
    113 {
    114   NOTIFY(TrajectoryChanged);
    115   AtomicPosition.pop_back();
    116   AtomicVelocity.pop_back();
    117   AtomicForce.pop_back();
     104void 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);
    118111  LOG(5,"AtomInfo::removeTrajectoryStep() called, size is ("
    119112      << AtomicPosition.size() << ","
    … …  
    141134const double& AtomInfo::operator[](size_t i) const
    142135{
    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];
     136  return atStep(i, WorldTime::getTime());
    148137}
    149138
    150139const double& AtomInfo::at(size_t i) const
    151140{
    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);
     141  return atStep(i, WorldTime::getTime());
    157142}
    158143
    159144const double& AtomInfo::atStep(size_t i, unsigned int _step) const
    160145{
    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);
     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];
    166151}
    167152
    … …  
    170155  OBSERVE;
    171156  NOTIFY(AtomObservable::PositionChanged);
    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;
     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 */
     179bool AtomInfo::isStepPresent(const unsigned int _step) const
     180{
     181  VectorTrajectory_t::const_iterator iter =
     182      AtomicPosition.find(_step);
     183  return iter != AtomicPosition.end();
    177184}
    178185
    179186const Vector& AtomInfo::getPosition() const
    180187{
    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()];
     188  return getPositionAtStep(WorldTime::getTime());
    186189}
    187190
    188191const Vector& AtomInfo::getPositionAtStep(const unsigned int _step) const
    189192{
    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];
     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;
    195198}
    196199
    … …  
    207210{
    208211  const element *elem = World::getInstance().getPeriode()->FindElement(Z);
    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 //}
     212  setType(elem);
     213}
    227214
    228215const Vector& AtomInfo::getAtomicVelocity() const
    229216{
    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()];
     217  return getAtomicVelocityAtStep(WorldTime::getTime());
    235218}
    236219
    237220const Vector& AtomInfo::getAtomicVelocityAtStep(const unsigned int _step) const
    238221{
    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];
     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;
    244231}
    245232
    246233void AtomInfo::setAtomicVelocity(const Vector &_newvelocity)
    247234{
     235  setAtomicVelocityAtStep(WorldTime::getTime(), _newvelocity);
     236}
     237
     238void AtomInfo::setAtomicVelocityAtStep(const unsigned int _step, const Vector &_newvelocity)
     239{
    248240  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;
    255 }
    256 
    257 void AtomInfo::setAtomicVelocityAtStep(const unsigned int _step, const Vector &_newvelocity)
    258 {
    259   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  }
    260253  if (WorldTime::getTime() == _step)
    261254    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   }
    273255}
    274256
    275257const Vector& AtomInfo::getAtomicForce() const
    276258{
    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()];
     259  return getAtomicForceAtStep(WorldTime::getTime());
    282260}
    283261
    284262const Vector& AtomInfo::getAtomicForceAtStep(const unsigned int _step) const
    285263{
    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];
     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;
    291273}
    292274
    293275void AtomInfo::setAtomicForce(const Vector &_newforce)
    294276{
     277  setAtomicForceAtStep(WorldTime::getTime(), _newforce);
     278}
     279
     280void AtomInfo::setAtomicForceAtStep(const unsigned int _step, const Vector &_newforce)
     281{
    295282  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;
    302 }
    303 
    304 void AtomInfo::setAtomicForceAtStep(const unsigned int _step, const Vector &_newforce)
    305 {
    306   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  }
    307295  if (WorldTime::getTime() == _step)
    308296    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   }
    320297}
    321298
    … …  
    334311void AtomInfo::setPosition(const Vector& _vector)
    335312{
     313  setPositionAtStep(WorldTime::getTime(), _vector);
     314}
     315
     316void AtomInfo::setPositionAtStep(unsigned int _step, const Vector& _vector)
     317{
    336318  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;
    344 }
    345 
    346 void AtomInfo::setPositionAtStep(unsigned int _step, const Vector& _vector)
    347 {
    348   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  }
    349331  if (WorldTime::getTime() == _step)
    350332    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;
    363333}
    364334
    365335const VectorInterface& AtomInfo::operator+=(const Vector& b)
    366336{
    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;
     337  setPosition(getPosition()+b);
    374338  return *this;
    375339}
    … …  
    377341const VectorInterface& AtomInfo::operator-=(const Vector& b)
    378342{
    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;
     343  setPosition(getPosition()-b);
    386344  return *this;
    387345}
    … …  
    389347Vector const AtomInfo::operator+(const Vector& b) const
    390348{
    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()]);
     349  Vector a(getPosition());
    396350  a += b;
    397351  return a;
    … …  
    400354Vector const AtomInfo::operator-(const Vector& b) const
    401355{
    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()]);
     356  Vector a(getPosition());
    407357  a -= b;
    408358  return a;
    … …  
    411361double AtomInfo::distance(const Vector &point) const
    412362{
    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);
     363  return getPosition().distance(point);
    418364}
    419365
    420366double AtomInfo::DistanceSquared(const Vector &y) const
    421367{
    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);
     368  return getPosition().DistanceSquared(y);
    427369}
    428370
    429371double AtomInfo::distance(const VectorInterface &_atom) const
    430372{
    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()]);
     373  return _atom.distance(getPosition());
    436374}
    437375
    438376double AtomInfo::DistanceSquared(const VectorInterface &_atom) const
    439377{
    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()]);
     378  return _atom.DistanceSquared(getPosition());
    445379}
    446380
    447381VectorInterface &AtomInfo::operator=(const Vector& _vector)
    448382{
    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;
     383  setPosition(_vector);
    456384  return *this;
    457385}
    … …  
    459387void AtomInfo::ScaleAll(const double *factor)
    460388{
    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);
     389  Vector temp(getPosition());
     390  temp.ScaleAll(factor);
     391  setPosition(temp);
    468392}
    469393
    470394void AtomInfo::ScaleAll(const Vector &factor)
    471395{
    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);
     396  Vector temp(getPosition());
     397  temp.ScaleAll(factor);
     398  setPosition(temp);
    479399}
    480400
    481401void AtomInfo::Scale(const double factor)
    482402{
    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);
     403  Vector temp(getPosition());
     404  temp.Scale(factor);
     405  setPosition(temp);
    490406}
    491407
    492408void AtomInfo::Zero()
    493409{
    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();
     410  setPosition(zeroVec);
    501411}
    502412
    503413void AtomInfo::One(const double one)
    504414{
    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);
     415  setPosition(Vector(one,one,one));
    512416}
    513417
    514418void AtomInfo::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
    515419{
    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);
     420  Vector newPos;
     421  newPos.LinearCombinationOfVectors(x1,x2,x3,factors);
     422  setPosition(newPos);
    523423}
    524424
    … …  
    528428double AtomInfo::getKineticEnergy(const unsigned int _step) const
    529429{
    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();
     430  return getMass() * getAtomicVelocityAtStep(_step).NormSquared();
    535431}
    536432
    537433Vector AtomInfo::getMomentum(const unsigned int _step) const
    538434{
    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.
     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 *
    548441 * \param MaxSteps
    549442 */
    550443void AtomInfo::ResizeTrajectory(size_t MaxSteps)
    551444{
    552   for (;AtomicPosition.size() <= (unsigned int)(MaxSteps);)
    553     UpdateSteps();
     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  }
    554464}
    555465
    556466size_t AtomInfo::getTrajectorySize() const
    557467{
    558   return AtomicPosition.size();
     468  // mind greater comp for map here: first element is latest in time steps!
     469  return AtomicPosition.begin()->first+1;
    559470}
    560471
    … …  
    562473{
    563474  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 */
     484void 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;
    564493}
    565494
    … …  
    579508  }
    580509
    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   }
     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);
    594525};
    595526
    … …  
    656587};
    657588
    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 
    678589std::ostream & AtomInfo::operator << (std::ostream &ost) const
    679590{

  • src/Atom/atom_atominfo.hpp

    rc7fe90 r5c8807  
    5656   * real functions, \sa AppendTrajectoryStep(), by all necessary subclasses.
    5757   */
    58   virtual void UpdateSteps()=0;
     58  virtual void UpdateStep(const unsigned int _step)=0;
    5959
    6060  /** Pops the last step in all trajectory vectors.
    … …  
    6464   * the real functions, \sa removeTrajectoryStep(), by all necessary subclasses.
    6565   */
    66   virtual void removeSteps()=0;
     66  virtual void removeStep(const unsigned int _step)=0;
    6767
    6868  /** DEPRECATED: Getter for element indicated by AtomicElement.
    … …  
    276276
    277277  // operations for trajectories
     278  bool isStepPresent(const unsigned int _step) const;
    278279  void ResizeTrajectory(size_t MaxSteps);
    279280  size_t getTrajectorySize() const;
    … …  
    297298   * vectors.
    298299   */
    299   void AppendTrajectoryStep();
     300  void AppendTrajectoryStep(const unsigned int _step);
    300301
    301302  /** Function used by this and inheriting classes to decrease the trajectory
    302303   * vectors by one.
    303304   */
    304   void removeTrajectoryStep();
     305  void removeTrajectoryStep(const unsigned int _step);
    305306
    306307  // make these protected only such that deriving atom class still has full
    307308  // access needed for clone and alike
    308   std::vector<Vector> AtomicPosition;       //!< coordinate vector of atom, giving last position within cell
    309   std::vector<Vector> AtomicVelocity;       //!< velocity vector of atom, giving last velocity within cell
    310   std::vector<Vector> AtomicForce;       //!< Force vector of atom, giving last force within cell
     309
     310  //!> typedef for a vector of Vectors with inverse sorting to make lower_bound return present or last past step
     311  typedef std::map<unsigned int, Vector, std::greater<unsigned int> > VectorTrajectory_t;
     312  VectorTrajectory_t AtomicPosition;       //!< coordinate vector of atom, giving last position within cell
     313  VectorTrajectory_t AtomicVelocity;       //!< velocity vector of atom, giving last velocity within cell
     314  VectorTrajectory_t AtomicForce;       //!< Force vector of atom, giving last force within cell
    311315
    312316private:

  • src/Atom/atom_bondedparticle.cpp

    rc7fe90 r5c8807  
    5252BondedParticle::BondedParticle()
    5353{
    54   ListOfBonds.push_back(BondList());
     54  ListOfBonds.insert( std::make_pair(0, emptyList) );
    5555};
    5656
    … …  
    5959BondedParticle::~BondedParticle()
    6060{
    61   removeAllBonds();
     61  for(BondTrajectory_t::iterator iter = ListOfBonds.begin(); !ListOfBonds.empty();
     62      iter = ListOfBonds.begin()) {
     63    removeAllBonds(iter->first);
     64    ListOfBonds.erase(iter);
     65  }
    6266};
    6367
    … …  
    220224}
    221225
    222 /** Removes all bonds in all timesteps and their instances, too.
     226/** Removes all bonds in current timestep and their instances, too.
    223227 *
    224228 */
    225229void BondedParticle::removeAllBonds()
    226230{
    227   for (size_t index = 0; index < ListOfBonds.size(); ++index)
    228     removeAllBonds(index);
     231  removeAllBonds(WorldTime::getTime());
    229232}
    230233
    … …  
    236239{
    237240  //LOG(3,"INFO: Clearing all bonds of " << *this << ": " << ListOfBonds[_step]);
    238   for (BondList::iterator iter = (ListOfBonds[_step]).begin();
    239       !(ListOfBonds[_step]).empty();
    240       iter = (ListOfBonds[_step]).begin()) {
    241     //LOG(3,"INFO: Clearing bond (" << *iter << ") " << *(*iter) << " of list " << &ListOfBonds);
    242     atom * const Other = (*iter)->GetOtherAtom(this);
    243     ASSERT( Other != NULL,
    244         "BondedParticle::removeAllBonds() - cannot find bond partner for "
    245         +toString(**iter)+".");
    246     Other->UnregisterBond(_step, *iter);
    247     UnregisterBond(_step, *iter);
    248   }
     241  BondTrajectory_t::iterator listiter = ListOfBonds.find(_step);
     242  if (listiter != ListOfBonds.end())
     243    for (BondList::iterator iter = listiter->second.begin();
     244        !listiter->second.empty();
     245        iter = listiter->second.begin()) {
     246      //LOG(3,"INFO: Clearing bond (" << *iter << ") " << *(*iter) << " of list " << &ListOfBonds);
     247      atom * const Other = (*iter)->GetOtherAtom(this);
     248      ASSERT( Other != NULL,
     249          "BondedParticle::removeAllBonds() - cannot find bond partner for "
     250          +toString(**iter)+".");
     251      Other->UnregisterBond(_step, *iter);
     252      UnregisterBond(_step, *iter);
     253    }
    249254}
    250255
    … …  
    260265    if (Binder->Contains(this)) {
    261266      //LOG(3,"INFO: Registering bond "<< *Binder << " with atom " << *this << " at step " << _step);
    262       if (ListOfBonds.size() <= _step)
    263         ListOfBonds.resize(_step+1);
    264       ListOfBonds[_step].push_back(Binder);
     267      std::pair< BondTrajectory_t::iterator, bool> inserter =
     268          ListOfBonds.insert( std::make_pair( _step, BondList(1, Binder)) );
     269      if (!inserter.second)
     270        inserter.first->second.push_back(Binder);
    265271      if (WorldTime::getTime() == _step)
    266272        NOTIFY(AtomObservable::BondsAdded);
    … …  
    290296      OBSERVE;
    291297      //LOG(0,"INFO: Unregistering bond "<< *Binder << " from list " << &ListOfBonds << " of atom " << *this << " at step " << step);
     298      BondTrajectory_t::iterator listiter = ListOfBonds.find(_step);
     299      if (listiter != ListOfBonds.end()) {
    292300#ifndef NDEBUG
    293       BondList::const_iterator iter =
    294           std::find(ListOfBonds[_step].begin(), ListOfBonds[_step].end(), Binder);
    295       ASSERT( iter != ListOfBonds[_step].end(),
    296           "BondedParticle::UnregisterBond() - "+toString(*Binder)+" not contained at "
    297           +toString(_step));
     301        BondList::const_iterator iter =
     302            std::find(listiter->second.begin(), listiter->second.end(), Binder);
     303        ASSERT( iter != listiter->second.end(),
     304            "BondedParticle::UnregisterBond() - "+toString(*Binder)+" not contained at "
     305            +toString(_step));
    298306#endif
    299       Binder->removeAtom(this);
    300       ListOfBonds[_step].remove(Binder);
    301       if (WorldTime::getTime() == _step)
    302         NOTIFY(AtomObservable::BondsRemoved);
    303       status = true;
     307        Binder->removeAtom(this);
     308        listiter->second.remove(Binder);
     309        if (WorldTime::getTime() == _step)
     310          NOTIFY(AtomObservable::BondsRemoved);
     311        status = true;
     312      }
    304313    } else {
    305314      ELOG(1, *Binder << " does not contain " << *this << ".");
    … …  
    328337{
    329338  int step = -1;
    330   int tempstep = 0;
    331   for(std::vector<BondList>::const_iterator iter = ListOfBonds.begin();
    332       iter != ListOfBonds.end();
    333       ++iter,++tempstep) {
    334     for (BondList::const_iterator bonditer = iter->begin();
    335         bonditer != iter->end();
     339  for(BondTrajectory_t::const_iterator listiter = ListOfBonds.begin();
     340      listiter != ListOfBonds.end();
     341      ++listiter) {
     342    for (BondList::const_iterator bonditer = listiter->second.begin();
     343        bonditer != listiter->second.end();
    336344        ++bonditer) {
    337345      if ((*bonditer) == Binder) {
    338         step = tempstep;
     346        step = listiter->first;
    339347        break;
    340348      }

  • src/Atom/atom_bondedparticleinfo.cpp

    rc7fe90 r5c8807  
    5757{}
    5858
    59 void BondedParticleInfo::AppendTrajectoryStep()
     59void BondedParticleInfo::AppendTrajectoryStep(const unsigned int _step)
    6060{
    61   ListOfBonds.push_back(BondList());
     61  ListOfBonds.insert( std::make_pair(_step, emptyList) );
    6262  LOG(5,"BondedParticleInfo::AppendTrajectoryStep() called, size is " << ListOfBonds.size());
    6363}
    6464
    65 void BondedParticleInfo::removeTrajectoryStep()
     65void BondedParticleInfo::removeTrajectoryStep(const unsigned int _step)
    6666{
    67   ListOfBonds.pop_back();
     67  ListOfBonds.erase(_step);
    6868  LOG(5,"BondedParticleInfo::removeTrajectoryStep() called, size is " << ListOfBonds.size());
    6969}
    … …  
    7171const BondList& BondedParticleInfo::getListOfBonds() const
    7272{
    73   if(WorldTime::getTime() < ListOfBonds.size())
    74     return ListOfBonds[WorldTime::getTime()];
    75   else
    76     return emptyList;
     73  return getListOfBondsAtStep(WorldTime::getTime());
    7774}
    78 
    79 //BondList& BondedParticleInfo::getListOfBonds()
    80 //{
    81 //  // todo: here we actually need a container on whose destruction notifiy is emitted, i.e.
    82 //  // similar or simply an ObservedContainer.
    83 //  OBSERVE;
    84 //  NOTIFY(AtomObservable::BondsAdded);
    85 //  const unsigned int size = ListOfBonds.size();
    86 //  ASSERT(WorldTime::getTime() <= size,
    87 //      "BondedParticleInfo::getBondsAtStep() - Access out of range: "
    88 //      +toString(WorldTime::getTime())
    89 //      +" not in [0,"+toString(size)+"].");
    90 //  if (WorldTime::getTime() == size) {
    91 //    UpdateSteps();
    92 //  }
    93 //  return ListOfBonds[WorldTime::getTime()];
    94 //}
    9575
    9676const BondList& BondedParticleInfo::getListOfBondsAtStep(unsigned int _step) const
    9777{
    98   if(_step < ListOfBonds.size())
    99     return ListOfBonds[_step];
    100   else
    101     return emptyList;
     78  BondTrajectory_t::const_iterator iter =
     79    ListOfBonds.find(_step);
     80  if (iter != ListOfBonds.end())
     81    return iter->second;
     82  return emptyList;
    10283}
    103 
    104 //BondList& BondedParticleInfo::getListOfBondsAtStep(unsigned int _step)
    105 //{
    106 //  const unsigned int size = ListOfBonds.size();
    107 //  ASSERT(_step <= size,
    108 //      "BondedParticleInfo::getBondsAtStep() - Access out of range: "
    109 //      +toString(_step)
    110 //      +" not in [0,"+toString(size)+"].");
    111 //  if (_step == size) {
    112 //    UpdateSteps();
    113 //  }
    114 //  return ListOfBonds[_step];
    115 //}

  • src/Atom/atom_bondedparticleinfo.hpp

    rc7fe90 r5c8807  
    4949   * real functions, \sa AppendTrajectoryStep(), by all necessary subclasses.
    5050   */
    51   virtual void UpdateSteps()=0;
     51  virtual void UpdateStep(const unsigned int _step)=0;
     52
     53  /** Pops the last step in all trajectory vectors.
     54   *
     55   * This allows to decrease all trajectories contained in different classes
     56   * by one consistently. This is implemented by the topmost class which calls
     57   * the real functions, \sa removeTrajectoryStep(), by all necessary subclasses.
     58   */
     59  virtual void removeStep(const unsigned int _step)=0;
    5260
    5361  /** Const accessor to ListOfBonds of WorldTime::CurrentTime.
    … …  
    8795   * vector.
    8896   */
    89   void AppendTrajectoryStep();
     97  void AppendTrajectoryStep(const unsigned int _step);
    9098
    9199  /** Function used by this and inheriting classes to reduce the ListOfBonds
    92100   * vector by one.
    93101   */
    94   void removeTrajectoryStep();
     102  void removeTrajectoryStep(const unsigned int _step);
    95103
    96   std::vector<BondList> ListOfBonds; //!< list of all bonds
     104  typedef std::map<unsigned int, BondList> BondTrajectory_t;
     105  BondTrajectory_t ListOfBonds; //!< list of all bonds
    97106  static BondList emptyList;  //!< empty list to return when step is not present
    98107};

  • src/Atom/unittests/AtomObserverUnitTest.cpp

    rc7fe90 r5c8807  
    130130void AtomObserverTest::AtomElementTest()
    131131{
    132   atoms[0]->setType(1);
     132  CPPUNIT_ASSERT( atoms[0]->getType()->getAtomicNumber() != 2 );
     133  atoms[0]->setType(2);
    133134
    134135  // check for update

  • src/Fragmentation/Exporters/HydrogenPool.cpp

    rc7fe90 r5c8807  
    103103  const size_t CurrentTime = WorldTime::getTime();
    104104  for (size_t step = _atom->getTrajectorySize(); step <= CurrentTime; ++step)
    105     _atom->UpdateSteps();
     105    _atom->UpdateStep(step);
    106106}
    107107

  • src/LinkedCell/unittests/stubs/TesselPointStub.cpp

    rc7fe90 r5c8807  
    4343  virtual ~TesselPoint();
    4444
    45   virtual void UpdateSteps();
     45  virtual void UpdateStep(const unsigned int _step);
    4646
    4747  TesselPoint * getTesselPoint();
    … …  
    5454{}
    5555
    56 void TesselPoint::UpdateSteps()
     56void TesselPoint::UpdateStep(const unsigned int _step)
    5757{}
    5858

  • src/Potentials/Makefile.am

    rc7fe90 r5c8807  
    5353  Potentials/Specifics/PotentialTypes.def \
    5454  Potentials/Specifics/PotentialTypes.undef
     55
     56if CONDLEVMAR
     57POTENTIALSHEADER += \
     58  Potentials/PotentialTrainer.hpp
     59POTENTIALSSOURCE = \
     60  Potentials/PotentialTrainer.cpp
     61endif
    5562
    5663# add here headers for real potentials (i.e. the PotentialFactory should instantiate) 

  • src/World.cpp

    rc7fe90 r5c8807  
    234234  if (_step != WorldTime::getTime()) {
    235235    const unsigned int oldstep = WorldTime::getTime();
    236     // set new time
     236
     237    // 1. copy bond graph (such not each addBond causes GUI update)
     238    if (!areBondsPresent(_step)) {
     239//      AtomComposite Set = getAllAtoms();
     240//      BG->cleanAdjacencyList(Set);
     241      copyBondgraph(oldstep, _step);
     242    }
     243
     244    // 2. set new time
    237245    WorldTime::getInstance().setTime(_step);
     246
    238247    // TODO: removed when BondGraph creates the adjacency
    239     // 1. remove all of World's molecules
     248    // 3. remove all of World's molecules
    240249    for (MoleculeIterator iter = getMoleculeIter();
    241250        getMoleculeIter() != moleculeEnd();
    … …  
    244253      destroyMolecule(*iter);
    245254    }
    246     // 2. copy bond graph
    247     if (!areBondsPresent(_step)) {
    248 //      AtomComposite Set = getAllAtoms();
    249 //      BG->cleanAdjacencyList(Set);
    250       copyBondgraph(oldstep, _step);
    251     }
    252 
    253     // 3. scan for connected subgraphs => molecules
     255
     256    // 4. scan for connected subgraphs => molecules
    254257    DepthFirstSearchAnalysis DFS;
    255258    DFS();

  • src/molecule_geometry.cpp

    rc7fe90 r5c8807  
    295295void molecule::Scale(const double *factor)
    296296{
    297   for (iterator iter = begin(); iter != end(); ++iter) {
    298     for (size_t j=0;j<(*iter)->getTrajectorySize();j++) {
    299       Vector temp = (*iter)->getPositionAtStep(j);
    300       temp.ScaleAll(factor);
    301       (*iter)->setPositionAtStep(j,temp);
    302     }
    303   }
     297  for (iterator iter = begin(); iter != end(); ++iter)
     298    for (size_t j=0;j<(*iter)->getTrajectorySize();j++)
     299      if ((*iter)->isStepPresent(j)) {
     300        Vector temp = (*iter)->getPositionAtStep(j);
     301        temp.ScaleAll(factor);
     302        (*iter)->setPositionAtStep(j,temp);
     303      }
    304304};
    305305
    … …  
    309309void molecule::Translate(const Vector &trans)
    310310{
    311   getAtomSet().translate(trans);
     311  for (iterator iter = begin(); iter != end(); ++iter)
     312    for (size_t j=0;j<(*iter)->getTrajectorySize();j++)
     313      if ((*iter)->isStepPresent(j))
     314        (*iter)->setPositionAtStep(j, (*iter)->getPositionAtStep(j) + (trans));
    312315};
    313316

  • tests/Python/AllActions/options.dat

    rc7fe90 r5c8807  
    9999mesh-size       "10,10,10"
    100100min-distance    "1."
     101mirror-atoms    "1.,1.,1."
    101102molecule-by-id  "0"
    102103near-field-cells        "3"
    … …  
    124125output-types    "xyz mpqc"
    125126parse-homologies        "homology.dat"
     127parse-potentials        "water.potentials"
    126128parse-tremolo-potentials        "argon.potentials"
    127129parse-tremolo-potentials        "tensid.potentials"
    … …  
    129131parser-parameters       "psi3"
    130132periodic        "0"
     133plane-offset    "5."
    131134position        "0 0 0"
    132135position        "0 0 1"
    … …  
    161164save-bonds      "test.bond"
    162165save-homologies "homology.dat"
     166save-potentials "water.potentials"
    163167save-selected-atoms     "testsave.xyz"
    164168save-selected-atoms-as-exttypes "test.exttypes"

  • tests/regression/Atoms/testsuite-atoms.at

    rc7fe90 r5c8807  
    4141# translate to origin
    4242m4_include([Atoms/TranslationToOrigin/testsuite-atoms-translation-to-origin.at])
     43
     44# mirror atoms
     45m4_include([Atoms/Mirror/testsuite-atoms-mirror.at])

  • tests/regression/Makefile.am

    rc7fe90 r5c8807  
    3535        $(srcdir)/Atoms/Add/testsuite-atoms-add.at \
    3636        $(srcdir)/Atoms/ChangeElement/testsuite-atoms-change-element.at \
     37        $(srcdir)/Atoms/Mirror/testsuite-atoms-mirror.at \
    3738        $(srcdir)/Atoms/Remove/testsuite-atoms-remove.at \
    3839        $(srcdir)/Atoms/RemoveCuboid/testsuite-atoms-remove-cuboid.at \

  • tests/regression/Potential/FitPotential/testsuite-potential-fit-potential.at

    rc7fe90 r5c8807  
    2626AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Potential/FitPotential/pre/$file $file], 0)
    2727AT_CHECK([chmod u+w $file], 0, [ignore], [ignore])
    28 AT_CHECK([../../molecuilder --parse-homologies $file --set-random-number-engine "lagged_fibonacci607" --random-number-engine-parameters "seed=1;" --set-random-number-distribution "uniform_real" --random-number-distribution-parameters "min=0;max=1;" --fit-potential --potential-type "morse" --potential-charges 8 1 --fragment-charges 1 8 1 --set-threshold 1e-6], 0, [stdout], [ignore])
     28AT_CHECK([../../molecuilder \
     29        --parse-homologies $file \
     30        --set-random-number-engine "lagged_fibonacci607" \
     31                --random-number-engine-parameters "seed=1;" \
     32        --set-random-number-distribution "uniform_real" \
     33                --random-number-distribution-parameters "min=0;max=1;" \
     34        --fit-potential \
     35                --potential-type "morse" \
     36                --potential-charges 8 1 \
     37                --fragment-charges 1 8 1 \
     38                --set-threshold 1e-6 \
     39        --save-potentials length.potentials], 0, [stdout], [ignore])
    2940# check that L_2 error is below 1e-6
    3041AT_CHECK([grep "Best parameters with L2 error" stdout | awk '{if ($8 > 1e-6) exit 1}'], 0, [ignore], [ignore])
    3142# check parameters to printed precision
    32 AT_CHECK([grep "morse:.*particle_type1=8,.*particle_type2=1,.*spring_constant=1.27.*,.*equilibrium_distance=1.78.*,.*dissociation_energy=0.19.*;" stdout], 0, [ignore], [ignore])
     43AT_CHECK([grep "morse:.*particle_type1=8,.*particle_type2=1,.*spring_constant=1.27.*,.*equilibrium_distance=1.78.*,.*dissociation_energy=0.19.*;" length.potentials], 0, [ignore], [ignore])
    3344
    3445AT_CLEANUP
    … …  
    4152AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Potential/FitPotential/pre/$file $file], 0)
    4253AT_CHECK([chmod u+w $file], 0, [ignore], [ignore])
    43 AT_CHECK([../../molecuilder --parse-homologies $file --set-random-number-engine "lagged_fibonacci607" --random-number-engine-parameters "seed=1;" --set-random-number-distribution "uniform_real" --random-number-distribution-parameters "min=0;max=1;" --fit-potential --potential-type "harmonic_bond" --potential-charges 8 1 --fragment-charges 1 8 1 --set-threshold 1e-6], 0, [stdout], [ignore])
     54AT_CHECK([../../molecuilder \
     55        --parse-homologies $file \
     56        --set-random-number-engine "lagged_fibonacci607" \
     57                --random-number-engine-parameters "seed=1;" \
     58        --set-random-number-distribution "uniform_real" \
     59                --random-number-distribution-parameters "min=0;max=1;" \
     60        --fit-potential \
     61                --potential-type "harmonic_bond" \
     62                --potential-charges 8 1 \
     63                --fragment-charges 1 8 1 \
     64                --set-threshold 1e-6 \
     65        --save-potentials harmonic.potentials], 0, [stdout], [ignore])
    4466# check that L_2 error is below 1e-6
    4567AT_CHECK([grep "Best parameters with L2 error" stdout | awk '{if ($8 > 1e-6) exit 1}'], 0, [ignore], [ignore])
    4668# check parameters to printed precision
    47 AT_CHECK([grep "harmonic_bond:.*particle_type1=8,.*particle_type2=1,.*spring_constant=0.29.*,.*equilibrium_distance=1.8.*;" stdout], 0, [ignore], [ignore])
     69AT_CHECK([grep "harmonic_bond:.*particle_type1=8,.*particle_type2=1,.*spring_constant=0.29.*,.*equilibrium_distance=1.8.*;" harmonic.potentials], 0, [ignore], [ignore])
    4870
    4971AT_CLEANUP
    … …  
    5678AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Potential/FitPotential/pre/$file $file], 0)
    5779AT_CHECK([chmod u+w $file], 0, [ignore], [ignore])
    58 AT_CHECK([../../molecuilder --parse-homologies $file --set-random-number-engine "lagged_fibonacci607" --random-number-engine-parameters "seed=1;" --set-random-number-distribution "uniform_real" --random-number-distribution-parameters "min=0;max=1;" --fit-potential --potential-type "harmonic_angle" --potential-charges 1 8 1 --fragment-charges 1 8 1 --set-threshold 1e-6], 0, [stdout], [ignore])
     80AT_CHECK([../../molecuilder \
     81        --parse-homologies $file \
     82        --set-random-number-engine "lagged_fibonacci607" \
     83                --random-number-engine-parameters "seed=1;" \
     84        --set-random-number-distribution "uniform_real" \
     85                --random-number-distribution-parameters "min=0;max=1;" \
     86        --fit-potential \
     87                --potential-type "harmonic_angle" \
     88                --potential-charges 1 8 1 \
     89                --fragment-charges 1 8 1 \
     90                --set-threshold 1e-6 \
     91        --save-potentials angle.potentials], 0, [stdout], [ignore])
    5992# check that L_2 error is below 1e-6
    6093AT_CHECK([grep "||e||_2:" stdout | awk '{if ($7 > 1e-6) exit 1}'], 0, [ignore], [ignore])
    6194# check parameters to printed precision
    62 AT_CHECK([grep "harmonic_angle:.*particle_type1=1,.*particle_type2=8,.*particle_type3=1,.*spring_constant=0.10.*,.*equilibrium_distance=-0.27.*;" stdout], 0, [ignore], [ignore])
     95AT_CHECK([grep "harmonic_angle:.*particle_type1=1,.*particle_type2=8,.*particle_type3=1,.*spring_constant=0.10.*,.*equilibrium_distance=-0.27.*;" angle.potentials], 0, [ignore], [ignore])
    6396
    6497AT_CLEANUP
    … …  
    71104AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Potential/FitPotential/pre/$file $file], 0)
    72105AT_CHECK([chmod u+w $file], 0, [ignore], [ignore])
    73 AT_CHECK([../../molecuilder --parse-homologies $file --set-random-number-engine "lagged_fibonacci607" --random-number-engine-parameters "seed=1;" --set-random-number-distribution "uniform_real" --random-number-distribution-parameters "min=0;max=1;" --fit-potential --potential-type "torsion" --potential-charges 6 6 6 6 --fragment-charges 6 6 6 6 1 1 1 1 1 1 1 1 1 1 --set-threshold 2e-10], 0, [stdout], [ignore])
     106AT_CHECK([../../molecuilder \
     107        --parse-homologies $file \
     108        --set-random-number-engine "lagged_fibonacci607" \
     109                --random-number-engine-parameters "seed=1;" \
     110        --set-random-number-distribution "uniform_real" \
     111                --random-number-distribution-parameters "min=0;max=1;" \
     112        --fit-potential \
     113                --potential-type "torsion" \
     114                --potential-charges 6 6 6 6 \
     115                --fragment-charges 6 6 6 6 1 1 1 1 1 1 1 1 1 1 \
     116                --set-threshold 2e-10 \
     117        --save-potentials torsion.potentials], 0, [stdout], [ignore])
    74118# check that L_2 error is below 9e-12 ... just 2e-10 otherwise test takes tooo long
    75119AT_CHECK([grep "Best parameters with L2 error" stdout | awk '{if ($8 > 2e-10) exit 1}'], 0, [ignore], [ignore])
    76 AT_CHECK([grep "torsion:.*particle_type1=6,.*particle_type2=6,.*particle_type3=6,.*particle_type4=6,.*spring_constant=.*,.*equilibrium_distance=.*;" stdout], 0, [ignore], [ignore])
    77 #AT_CHECK([grep "torsion:.*particle_type1=6,.*particle_type2=6,.*particle_type3=6,.*particle_type4=6,.*spring_constant=0.001.*,.*equilibrium_distance=0.99.*;" stdout], 0, [ignore], [ignore])
     120AT_CHECK([grep "torsion:.*particle_type1=6,.*particle_type2=6,.*particle_type3=6,.*particle_type4=6,.*spring_constant=.*,.*equilibrium_distance=.*;" torsion.potentials], 0, [ignore], [ignore])
     121#AT_CHECK([grep "torsion:.*particle_type1=6,.*particle_type2=6,.*particle_type3=6,.*particle_type4=6,.*spring_constant=0.001.*,.*equilibrium_distance=0.99.*;" torsion.potentials], 0, [ignore], [ignore])
    78122
    79123AT_CLEANUP
    … …  
    86130AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Potential/FitPotential/pre/$file $file], 0)
    87131AT_CHECK([chmod u+w $file], 0, [ignore], [ignore])
    88 AT_CHECK([../../molecuilder --parse-homologies $file --set-random-number-engine "lagged_fibonacci607" --random-number-engine-parameters "seed=1;" --set-random-number-distribution "uniform_real" --random-number-distribution-parameters "min=0;max=1;" --fit-potential --potential-type "improper" --potential-charges 1 7 1 1 --fragment-charges 7 1 1 1 --set-threshold 3e-4], 0, [stdout], [ignore])
     132AT_CHECK([../../molecuilder \
     133        --parse-homologies $file \
     134        --set-random-number-engine "lagged_fibonacci607" \
     135                --random-number-engine-parameters "seed=1;" \
     136        --set-random-number-distribution "uniform_real" \
     137                --random-number-distribution-parameters "min=0;max=1;" \
     138        --fit-potential \
     139                --potential-type "improper" \
     140                --potential-charges 1 7 1 1 \
     141                --fragment-charges 7 1 1 1 \
     142                --set-threshold 3e-4 \
     143        --save-potentials improper.potentials], 0, [stdout], [ignore])
    89144# check that L_2 error is below 3e-4
    90145AT_CHECK([grep "Best parameters with L2 error" stdout | awk '{if ($8 > 3e-4) exit 1}'], 0, [ignore], [ignore])
    91146# check parameters to printed precision
    92 AT_CHECK([grep "improper:.*particle_type1=1,.*particle_type2=7,.*particle_type3=1,.*particle_type4=1,.*spring_constant=.*,.*equilibrium_distance=.*;" stdout], 0, [ignore], [ignore])
    93 #AT_CHECK([grep "improper:.*particle_type1=1,.*particle_type2=7,.*particle_type3=1,.*particle_type4=1,.*spring_constant=1.02.*,.*equilibrium_distance=0.85.*;" stdout], 0, [ignore], [ignore])
     147AT_CHECK([grep "improper:.*particle_type1=1,.*particle_type2=7,.*particle_type3=1,.*particle_type4=1,.*spring_constant=.*,.*equilibrium_distance=.*;" improper.potentials], 0, [ignore], [ignore])
     148#AT_CHECK([grep "improper:.*particle_type1=1,.*particle_type2=7,.*particle_type3=1,.*particle_type4=1,.*spring_constant=1.02.*,.*equilibrium_distance=0.85.*;" improper.potentials], 0, [ignore], [ignore])
    94149
    95150AT_CLEANUP
    … …  
    102157AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Potential/FitPotential/pre/$file $file], 0)
    103158AT_CHECK([chmod u+w $file], 0, [ignore], [ignore])
    104 AT_CHECK([../../molecuilder --parse-homologies $file --set-random-number-engine "lagged_fibonacci607" --random-number-engine-parameters "seed=5;" --set-random-number-distribution "uniform_real" --random-number-distribution-parameters "min=0;max=1;" --fit-potential --potential-type "lennardjones" --potential-charges 18 18 --fragment-charges 18 18 --set-threshold 7e-9], 0, [stdout], [ignore])
     159AT_CHECK([../../molecuilder \
     160        --parse-homologies $file \
     161        --set-random-number-engine "lagged_fibonacci607" \
     162                --random-number-engine-parameters "seed=5;" \
     163        --set-random-number-distribution "uniform_real" \
     164                --random-number-distribution-parameters "min=0;max=1;" \
     165        --fit-potential \
     166                --potential-type "lennardjones" \
     167                --potential-charges 18 18 \
     168                --fragment-charges 18 18 \
     169                --set-threshold 7e-9 \
     170        --save-potentials lj.potentials], 0, [stdout], [ignore])
    105171# check that L_2 error is below 7e-11 ... just 7e-9 otherwise test takes too long
    106172AT_CHECK([grep "Best parameters with L2 error" stdout | awk '{if ($8 > 7e-9) exit 1}'], 0, [ignore], [ignore])
    107173# check parameters to printed precision
    108 AT_CHECK([grep "lennardjones:.*particle_type1=18,.*particle_type2=18,.*epsilon=.*,.*sigma=.*;" stdout], 0, [ignore], [ignore])
    109 #AT_CHECK([grep "lennardjones:.*particle_type1=18,.*particle_type2=18,.*epsilon=1.*e-05,.*sigma=8.2.*;" stdout], 0, [ignore], [ignore])
     174AT_CHECK([grep "lennardjones:.*particle_type1=18,.*particle_type2=18,.*epsilon=.*,.*sigma=.*;" lj.potentials], 0, [ignore], [ignore])
     175#AT_CHECK([grep "lennardjones:.*particle_type1=18,.*particle_type2=18,.*epsilon=1.*e-05,.*sigma=8.2.*;" lj.potentials], 0, [ignore], [ignore])
    110176
    111177AT_CLEANUP
Note: See TracChangeset for help on using the changeset viewer.