Changeset 83956e for src/Dynamics


Ignore:
Timestamp:
Jun 20, 2018, 8:20:42 AM (7 years ago)
Author:
Frederik Heber <frederik.heber@…>
Branches:
AutomationFragmentation_failures, Candidate_v1.6.1, ChemicalSpaceEvaluator, Exclude_Hydrogens_annealWithBondGraph, ForceAnnealing_with_BondGraph_contraction-expansion, StoppableMakroAction
Children:
646f73
Parents:
6a5921
git-author:
Frederik Heber <frederik.heber@…> (08/10/17 15:35:45)
git-committer:
Frederik Heber <frederik.heber@…> (06/20/18 08:20:42)
Message:

Rewrote annealWithBondGraph_BarzilaiBorwein() to simply distinguish expansion and contraction in bonds.

  • we shift the neighboring set away in case of expansion and towards in case of contraction.
Location:
src/Dynamics
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • src/Dynamics/BondVectors.cpp

    r6a5921 r83956e  
    285285}
    286286
     287void BondVectors::getProjectedGradientsForAtomAtStep(
     288    const atom &_walker,
     289    const Vector &_walkerGradient,
     290    const size_t _timestep,
     291    std::vector< std::vector<double> > &_projected_forces) const
     292{
     293  // gather subset of BondVectors for the current atom
     294  const std::vector<Vector> BondVectors = getAtomsBondVectorsAtStep(_walker, _timestep);
     295  const BondList& ListOfBonds = _walker.getListOfBonds(); // we always use current bonds
     296
     297  // go over all its bonds
     298  std::vector<Vector>::const_iterator vectoriter = BondVectors.begin();
     299  for(BondList::const_iterator bonditer = ListOfBonds.begin();
     300      bonditer != ListOfBonds.end(); ++bonditer, ++vectoriter) {
     301    const bond::ptr &current_bond = *bonditer;
     302    const Vector &BondVector = *vectoriter;
     303
     304    // bv goes from rightatom to leftatom
     305    // as plus sign in force indicates expansion, minus indicates contraction
     306    // leftatom: from right to left means same sign expansion, opposite contraction
     307    // rightatom: from right to left means opposite sign expansion, same sign contraction
     308    const double sign = (&_walker == current_bond->leftatom) ? 1. : -1.;
     309    const double temp = sign*_walkerGradient.ScalarProduct(BondVector);
     310    LOG(4, "DEBUG: BondVector " << BondVector << " receives projected force of "
     311        << sign << "*" << _walkerGradient << "*" << BondVector << " = " << temp);
     312    const size_t index = getIndexForBond(current_bond);
     313    ASSERT( index != (size_t)-1,
     314        "ForceAnnealing() - could not find bond "+toString(*current_bond)
     315        +" in bondvectors");
     316    std::vector<double> &forcelist = (&_walker == current_bond->leftatom) ?
     317        _projected_forces[leftside] : _projected_forces[rightside];
     318    forcelist[index] = temp;
     319  }
     320  ASSERT( vectoriter == BondVectors.end(),
     321      "BondVectors::getRemnantGradientForAtomAtStep() - vectoriter is not at end when it should be.");
     322}
     323
    287324Vector BondVectors::getRemnantGradientForAtomAtStep(
    288325    const atom &_walker,
  • src/Dynamics/BondVectors.hpp

    r6a5921 r83956e  
    5656      const size_t &_step);
    5757
     58  /** Returns the number of bonds whose bond vectors have been registered.
     59   *
     60   * \param number of bonds
     61   */
     62  size_t size() const
     63  { return container.size(); }
     64
    5865  /** Getter for the sorted bonds.
    5966   *
     
    144151      forcestore_t _forcestore) const;
    145152
     153  /** Calculates the \a _walkkerGradient projected onto the bond vector for
     154   * every of \a _walker's bonds.
     155   *
     156   * \param _walker atom
     157   * \param _walkerGradient atom's gradient
     158   * \param _timestep time step
     159   * \param _projected_forces list of forces for every \a whichatom_t and every bond
     160   */
     161  void getProjectedGradientsForAtomAtStep(
     162      const atom &_walker,
     163      const Vector &_walkerGradient,
     164      const size_t _timestep,
     165      std::vector< std::vector<double> > &_projected_forces) const;
     166
     167  // knowing the number of bonds in total, we can setup the storage for the
     168  // projected forces
     169  enum whichatom_t {
     170    leftside=0,
     171    rightside=1,
     172    MAX_sides
     173  };
     174
    146175private:
    147176  /** Calculates the bond vector for each bond in the internal container.
     
    177206  //!> internal map for bond Bondvector association
    178207  mutable mapped_t current_mapped_vectors;
     208
    179209};
    180210
  • src/Dynamics/ForceAnnealing.hpp

    r6a5921 r83956e  
    207207            "ForceAnnealing::anneal() - if currentStep is "+toString(currentStep)
    208208            +", then there should be at least three time steps.");
    209         const Vector oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
     209        const Vector &oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
    210210        const Vector PositionDifference = currentPosition - oldPosition;
    211211        LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
     
    259259        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
    260260      // atom's force vector gives steepest descent direction
    261       const Vector oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
    262       const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
    263       const Vector oldGradient = (*iter)->getAtomicForceAtStep(OldTimeStep);
    264       const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
     261      const Vector &oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
     262      const Vector &currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
     263      const Vector &oldGradient = (*iter)->getAtomicForceAtStep(OldTimeStep);
     264      const Vector &currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
    265265      LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
    266266      LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition);
     
    304304  }
    305305
    306   // knowing the number of bonds in total, we can setup the storage for the
    307   // projected forces
    308   enum whichatom_t {
    309     leftside=0,
    310     rightside=1,
    311     MAX_sides
    312   };
    313 
    314   /** Helper function to put bond force into a container.
    315    *
    316    * \param _walker atom
    317    * \param _current_bond current bond of \a _walker
    318    * \param _timestep time step
    319    * \param _force calculated bond force
    320    * \param _bv bondvectors for obtaining the correct index
    321    * \param _projected_forces container
    322    */
    323   static void ForceStore(
    324       const atom &_walker,
    325       const bond::ptr &_current_bond,
    326       const size_t &_timestep,
    327       const double _force,
    328       const BondVectors &_bv,
    329       std::vector< // time step
    330             std::vector< // which bond side
    331               std::vector<double> > // over all bonds
    332                 > &_projected_forces)
    333   {
    334     std::vector<double> &forcelist = (&_walker == _current_bond->leftatom) ?
    335         _projected_forces[_timestep][leftside] : _projected_forces[_timestep][rightside];
    336     const size_t index = _bv.getIndexForBond(_current_bond);
    337     ASSERT( index != (size_t)-1,
    338         "ForceAnnealing() - could not find bond "+toString(*_current_bond)
    339         +" in bondvectors");
    340     forcelist[index] = _force;
    341   }
    342 
    343306  /** Performs Gradient optimization on the bonds with BarzilaiBorwein stepwdith.
    344307   *
     
    377340    BreadthFirstSearchGatherer NodeGatherer(BGcreator);
    378341
    379     /// We assume that a force is local, i.e. a bond is too short yet and hence
    380     /// the atom needs to be moved. However, all the adjacent (bound) atoms might
    381     /// already be at the perfect distance. If we just move the atom alone, we ruin
    382     /// all the other bonds. Hence, it would be sensible to move every atom found
    383     /// through the bond graph in the direction of the force as well by the same
    384     /// PositionUpdate. This is almost what we are going to do.
    385 
    386     /// One issue is: If we need to shorten bond, then we use the PositionUpdate
    387     /// also on the the other bond partner already. This is because it is in the
    388     /// direction of the bond. Therefore, the update is actually performed twice on
    389     /// each bond partner, i.e. the step size is twice as large as it should be.
    390     /// This problem only occurs when bonds need to be shortened, not when they
    391     /// need to be made longer (then the force vector is facing the other
    392     /// direction than the bond vector).
    393     /// As a remedy we need to average the force on either end of the bond and
    394     /// check whether each gradient points inwards out or outwards with respect
    395     /// to the bond and then shift accordingly.
    396 
    397     /// One more issue is that the projection onto the bond directions does not
    398     /// recover the gradient but may be larger as the bond directions are a
    399     /// generating system and not a basis (e.g. 3 bonds on a plane where 2 would
    400     /// suffice to span the plane). To this end, we need to account for the
    401     /// overestimation and obtain a weighting for each bond.
     342    /** We assume that a force is local, i.e. a bond is too short yet and hence
     343     * the atom needs to be moved. However, all the adjacent (bound) atoms might
     344     * already be at the perfect distance. If we just move the atom alone, we ruin
     345     * all the other bonds. Hence, it would be sensible to move every atom found
     346     * through the bond graph in the direction of the force as well by the same
     347     * PositionUpdate. This is almost what we are going to do, see below.
     348     *
     349     * This is to make the force a little more global in the sense of a multigrid
     350     * solver that uses various coarser grids to transport errors more effectively
     351     * over finely resolved grids.
     352     *
     353     */
     354
     355    /** The idea is that we project the gradients onto the bond vectors and determine
     356     * from the sum of projected gradients from either side whether the bond is
     357     * to contract or to expand. As the gradient acting as the normal vector of
     358     * a plane supported at the position of the atom separates all bonds into two
     359     * sets, we check whether all on one side are contracting and all on the other
     360     * side are expanding. In this case we may move not only the atom itself but
     361     * may propagate its update along a limited-horizon BFS to neighboring atoms.
     362     *
     363     */
    402364
    403365    // initialize helper class for bond vectors using bonds from range of atoms
     
    407369        AtomicForceManipulator<T>::atoms.end(),
    408370        _TimeStep); // use time step to update here as this is the current set of bonds
    409     const BondVectors::container_t &sorted_bonds = bv.getSorted();
    410 
    411     std::vector< // time step
    412       std::vector< // which bond side
    413         std::vector<double> > // over all bonds
    414           > projected_forces(2); // one for leftatoms, one for rightatoms (and for both time steps)
    415     for (size_t i=0;i<2;++i) {
    416       projected_forces[i].resize(MAX_sides);
    417       for (size_t j=0;j<MAX_sides;++j)
    418         projected_forces[i][j].resize(sorted_bonds.size(), 0.);
    419     }
    420 
    421     // for each atom we need to gather weights and then project the gradient
    422     typedef std::map<atomId_t, BondVectors::weights_t > weights_per_atom_t;
    423     std::vector<weights_per_atom_t> weights_per_atom(2);
    424     typedef std::map<atomId_t, Vector> RemnantGradient_per_atom_t;
    425     RemnantGradient_per_atom_t RemnantGradient_per_atom;
    426     for (size_t timestep = 0; timestep <= 1; ++timestep) {
    427       const size_t ReferenceTimeStep = CurrentTimeStep-timestep;
    428       LOG(2, "DEBUG: given time step is " << _TimeStep
    429           << ", timestep is " << timestep
    430           << ", and ReferenceTimeStep is " << ReferenceTimeStep);
    431 
    432       for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin();
    433           iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
    434         const atom &walker = *(*iter);
    435         const Vector &walkerGradient = walker.getAtomicForceAtStep(ReferenceTimeStep);
    436         LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
    437             << walker << " is " << walkerGradient << " with magnitude of "
    438             << walkerGradient.Norm());
    439 
    440         const BondList& ListOfBonds = walker.getListOfBonds();
    441         if (walkerGradient.Norm() > MYEPSILON) {
    442 
    443           // gather subset of BondVectors for the current atom
    444           const std::vector<Vector> BondVectors =
    445               bv.getAtomsBondVectorsAtStep(walker, ReferenceTimeStep);
    446 
    447           // go through all its bonds and calculate what magnitude is represented
    448           // by the others i.e. sum of scalar products against other bonds
    449           const std::pair<weights_per_atom_t::iterator, bool> inserter =
    450               weights_per_atom[timestep].insert(
    451                   std::make_pair(walker.getId(),
    452                   bv.getWeightsForAtomAtStep(walker, BondVectors, ReferenceTimeStep)) );
    453           ASSERT( inserter.second,
    454               "ForceAnnealing::operator() - weight map for atom "+toString(walker)
    455               +" and time step "+toString(timestep)+" already filled?");
    456           BondVectors::weights_t &weights = inserter.first->second;
    457           ASSERT( weights.size() == ListOfBonds.size(),
    458               "ForceAnnealing::operator() - number of weights "
    459               +toString(weights.size())+" does not match number of bonds "
    460               +toString(ListOfBonds.size())+", error in calculation?");
    461 
    462           // projected gradient over all bonds and place in one of projected_forces
    463           // using the obtained weights
    464           BondVectors::forcestore_t forcestoring =
    465               boost::bind(&ForceAnnealing::ForceStore, _1, _2, _3, _4,
    466                   boost::cref(bv), boost::ref(projected_forces));
    467           const Vector RemnantGradient = bv.getRemnantGradientForAtomAtStep(
    468               walker, walkerGradient, BondVectors, weights, timestep, forcestoring
    469           );
    470           RemnantGradient_per_atom.insert( std::make_pair(walker.getId(), RemnantGradient) );
    471         } else {
    472           LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
    473               << MYEPSILON << " for atom " << walker);
    474           // note that projected_forces is initialized to full length and filled
    475           // with zeros. Hence, nothing to do here
     371
     372    std::vector< // which bond side
     373      std::vector<double> > // over all bonds
     374        projected_forces; // one for leftatoms, one for rightatoms
     375    projected_forces.resize(BondVectors::MAX_sides);
     376    for (size_t j=0;j<BondVectors::MAX_sides;++j)
     377      projected_forces[j].resize(bv.size(), 0.);
     378
     379    // for each atom we need to project the gradient
     380    for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin();
     381        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
     382      const atom &walker = *(*iter);
     383      const Vector &walkerGradient = walker.getAtomicForceAtStep(CurrentTimeStep);
     384      const double GradientNorm = walkerGradient.Norm();
     385      LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
     386          << walker << " is " << walkerGradient << " with magnitude of "
     387          << GradientNorm);
     388
     389      if (GradientNorm > MYEPSILON) {
     390        bv.getProjectedGradientsForAtomAtStep(
     391            walker, walkerGradient, CurrentTimeStep, projected_forces
     392        );
     393      } else {
     394        LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
     395            << MYEPSILON << " for atom " << walker);
     396        // note that projected_forces is initialized to full length and filled
     397        // with zeros. Hence, nothing to do here
     398      }
     399    }
     400
     401    std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
     402    for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
     403        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
     404      atom &walker = *(*iter);
     405
     406      /// calculate step width
     407      const Vector &oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
     408      const Vector &currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
     409      const Vector &oldGradient = (*iter)->getAtomicForceAtStep(OldTimeStep);
     410      const Vector &currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
     411      LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
     412      LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition);
     413      LOG(4, "DEBUG: oldGradient for atom #" << (*iter)->getId() << " is " << oldGradient);
     414      LOG(4, "DEBUG: currentGradient for atom #" << (*iter)->getId() << " is " << currentGradient);
     415//      LOG(4, "DEBUG: Force for atom #" << (*iter)->getId() << " is " << currentGradient);
     416
     417      // we use Barzilai-Borwein update with position reversed to get descent
     418      const Vector PositionDifference = currentPosition - oldPosition;
     419      const Vector GradientDifference = (currentGradient - oldGradient);
     420      double stepwidth = 0.;
     421      if (GradientDifference.Norm() > MYEPSILON)
     422        stepwidth = fabs(PositionDifference.ScalarProduct(GradientDifference))/
     423            GradientDifference.NormSquared();
     424      if (fabs(stepwidth) < 1e-10) {
     425        // dont' warn in first step, deltat usage normal
     426        if (currentStep != 1)
     427          ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
     428        stepwidth = currentDeltat;
     429      }
     430      Vector PositionUpdate = stepwidth * currentGradient;
     431      LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
     432
     433      /** go through all bonds and check projected_forces and side of plane
     434       * the idea is that if all bonds on one side are contracting ones or expanding,
     435       * respectively, then we may shift not only the atom with respect to its
     436       * gradient but also its neighbors (towards contraction or towards
     437       * expansion depending on direction of gradient).
     438       * if they are mixed on both sides of the plane, then we simply shift
     439       * only the atom itself.
     440       * if they are not mixed on either side, then we also only shift the
     441       * atom, namely away from expanding and towards contracting bonds.
     442       */
     443
     444      // sign encodes side of plane and also encodes contracting(-) or expanding(+)
     445      typedef std::vector<int> sides_t;
     446      typedef std::vector<int> types_t;
     447      sides_t sides;
     448      types_t types;
     449      const BondList& ListOfBonds = walker.getListOfBonds();
     450      for(BondList::const_iterator bonditer = ListOfBonds.begin();
     451          bonditer != ListOfBonds.end(); ++bonditer) {
     452        const bond::ptr &current_bond = *bonditer;
     453
     454        // BondVector goes from bond::rightatom to bond::leftatom
     455        const size_t index = bv.getIndexForBond(current_bond);
     456        std::vector<double> &forcelist = (&walker == current_bond->leftatom) ?
     457            projected_forces[BondVectors::leftside] : projected_forces[BondVectors::rightside];
     458        const double &temp = forcelist[index];
     459        sides.push_back( -1.*temp/fabs(temp) ); // BondVectors has exactly opposite sign for sides decision
     460        const double sum =
     461            projected_forces[BondVectors::leftside][index]+projected_forces[BondVectors::rightside][index];
     462        types.push_back( sum/fabs(sum) );
     463        LOG(4, "DEBUG: Bond " << *current_bond << " is on side " << sides.back()
     464            << " and has type " << types.back());
     465      }
     466//      /// check whether both conditions are compatible:
     467//      // i.e. either we have ++/-- for all entries in sides and types
     468//      // or we have +-/-+ for all entries
     469//      // hence, multiplying and taking the sum and its absolute value
     470//      // should be equal to the maximum number of entries
     471//      sides_t results;
     472//      std::transform(
     473//          sides.begin(), sides.end(),
     474//          types.begin(),
     475//          std::back_inserter(results),
     476//          std::multiplies<int>);
     477//      int result = abs(std::accumulate(results.begin(), results.end(), 0, std::plus<int>));
     478
     479      std::vector<size_t> first_per_side(2, (size_t)-1); //!< mark down one representative from either side
     480      std::vector< std::vector<int> > types_per_side(2); //!< gather all types on each side
     481      types_t::const_iterator typesiter = types.begin();
     482      for (sides_t::const_iterator sidesiter = sides.begin();
     483          sidesiter != sides.end(); ++sidesiter, ++typesiter) {
     484        const size_t index = (*sidesiter+1)/2;
     485        types_per_side[index].push_back(*typesiter);
     486        if (first_per_side[index] == (size_t)-1)
     487          first_per_side[index] = std::distance(const_cast<const sides_t &>(sides).begin(), sidesiter);
     488      }
     489      LOG(4, "DEBUG: First on side minus is " << first_per_side[0] << ", and first on side plus is "
     490          << first_per_side[1]);
     491      //!> enumerate types per side with a little witching with the numbers to allow easy setting from types
     492      enum whichtypes_t {
     493        contracting=0,
     494        unset=1,
     495        expanding=2,
     496        mixed
     497      };
     498      std::vector<int> typeside(2, unset);
     499      for(size_t i=0;i<2;++i) {
     500        for (std::vector<int>::const_iterator tpsiter = types_per_side[i].begin();
     501            tpsiter != types_per_side[i].end(); ++tpsiter) {
     502          if (typeside[i] == unset) {
     503            typeside[i] = *tpsiter+1; //contracting(0) or expanding(2)
     504          } else {
     505            if (typeside[i] != (*tpsiter+1)) // no longer he same type
     506              typeside[i] = mixed;
     507          }
    476508        }
    477509      }
    478     }
    479 
    480     // step through each bond and shift the atoms
    481     std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
    482 
    483     LOG(3, "DEBUG: current step is " << currentStep << ", given time step is " << CurrentTimeStep);
    484     const BondVectors::mapped_t bondvectors = bv.getBondVectorsAtStep(CurrentTimeStep);
    485 
    486     for (BondVectors::container_t::const_iterator bondsiter = sorted_bonds.begin();
    487         bondsiter != sorted_bonds.end(); ++bondsiter) {
    488       const bond::ptr &current_bond = *bondsiter;
    489       const size_t index = bv.getIndexForBond(current_bond);
    490       const atom* bondatom[MAX_sides] = {
    491           current_bond->leftatom,
    492           current_bond->rightatom
    493       };
    494 
    495       // remove the edge
     510      LOG(4, "DEBUG: Minus side is " << typeside[0] << " and plus side is " << typeside[1]);
     511
     512      typedef std::vector< std::pair<atomId_t, atomId_t> > RemovedEdges_t;
     513      if ((typeside[0] != mixed) || (typeside[1] != mixed)) {
     514        const size_t sideno = ((typeside[0] != mixed) && (typeside[0] != unset)) ? 0 : 1;
     515        LOG(4, "DEBUG: Chosen side is " << sideno << " with type " << typeside[sideno]);
     516        ASSERT( (typeside[sideno] == contracting) || (typeside[sideno] == expanding),
     517            "annealWithBondGraph_BB() - chosen side is either expanding nor contracting.");
     518        // one side is not mixed, all bonds on one side are of same type
     519        // hence, find out which bonds to exclude
     520        const BondList& ListOfBonds = walker.getListOfBonds();
     521
     522        // away from gradient (minus) and contracting
     523        // or towards gradient (plus) and expanding
     524        // gather all on same side and remove
     525        const double sign =
     526            (sides[first_per_side[sideno]] == types[first_per_side[sideno]])
     527            ? sides[first_per_side[sideno]] : -1.*sides[first_per_side[sideno]];
     528        LOG(4, "DEBUG: Removing edges from side with sign " << sign);
     529        BondList::const_iterator bonditer = ListOfBonds.begin();
     530        RemovedEdges_t RemovedEdges;
     531        for (sides_t::const_iterator sidesiter = sides.begin();
     532            sidesiter != sides.end(); ++sidesiter, ++bonditer) {
     533          if (*sidesiter == sign) {
     534            // remove the edge
     535            const bond::ptr &current_bond = *bonditer;
     536            LOG(5, "DEBUG: Removing edge " << *current_bond);
     537            RemovedEdges.push_back( std::make_pair(
     538                current_bond->leftatom->getId(),
     539                current_bond->rightatom->getId())
     540            );
    496541#ifndef NDEBUG
    497       const bool status =
     542            const bool status =
    498543#endif
    499           BGcreator.removeEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
    500       ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
    501 
    502       // gather nodes for either atom
    503       BoostGraphHelpers::Nodeset_t bondside_set[MAX_sides];
    504       BreadthFirstSearchGatherer::distance_map_t distance_map[MAX_sides];
    505       for (size_t side=leftside;side<MAX_sides;++side) {
    506         bondside_set[side] = NodeGatherer(bondatom[side]->getId(), max_distance);
    507         distance_map[side] = NodeGatherer.getDistances();
    508         std::sort(bondside_set[side].begin(), bondside_set[side].end());
    509       }
    510 
    511       // re-add edge
    512       BGcreator.addEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
    513 
    514       // do for both leftatom and rightatom of bond
    515       for (size_t side = leftside; side < MAX_sides; ++side) {
    516         const double &bondforce = projected_forces[0][side][index];
    517         const double &oldbondforce = projected_forces[1][side][index];
    518         const double bondforcedifference = fabs(bondforce - oldbondforce);
    519         LOG(4, "DEBUG: bondforce for " << (side == leftside ? "left" : "right")
    520             << " side of bond is " << bondforce);
    521         LOG(4, "DEBUG: oldbondforce for " << (side == leftside ? "left" : "right")
    522             << " side of bond is " << oldbondforce);
    523         // if difference or bondforce itself is zero, do nothing
    524         if ((fabs(bondforce) < MYEPSILON) || (fabs(bondforcedifference) < MYEPSILON))
    525           continue;
    526 
    527         // get BondVector to bond
    528         const BondVectors::mapped_t::const_iterator bviter =
    529             bondvectors.find(current_bond);
    530         ASSERT( bviter != bondvectors.end(),
    531             "ForceAnnealing() - cannot find current_bond ?");
    532         ASSERT( fabs(bviter->second.Norm() -1.) < MYEPSILON,
    533             "ForceAnnealing() - norm of BondVector is not one");
    534         const Vector &BondVector = bviter->second;
    535 
    536         // calculate gradient and position differences for stepwidth
    537         const Vector currentGradient = bondforce * BondVector;
    538         LOG(4, "DEBUG: current projected gradient for "
    539             << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient);
    540         const Vector &oldPosition = bondatom[side]->getPositionAtStep(OldTimeStep);
    541         const Vector &currentPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep);
    542         const Vector PositionDifference = currentPosition - oldPosition;
    543         LOG(4, "DEBUG: old position is " << oldPosition);
    544         LOG(4, "DEBUG: current position is " << currentPosition);
    545         LOG(4, "DEBUG: difference in position is " << PositionDifference);
    546         LOG(4, "DEBUG: bondvector is " << BondVector);
    547         const double projected_PositionDifference = PositionDifference.ScalarProduct(BondVector);
    548         LOG(4, "DEBUG: difference in position projected onto bondvector is "
    549             << projected_PositionDifference);
    550         LOG(4, "DEBUG: abs. difference in forces is " << bondforcedifference);
    551 
    552         // calculate step width
    553         double stepwidth =
    554             fabs(projected_PositionDifference)/bondforcedifference;
    555         if (fabs(stepwidth) < 1e-10) {
    556           // dont' warn in first step, deltat usage normal
    557           if (currentStep != 1)
    558             ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
    559           stepwidth = currentDeltat;
     544                BGcreator.removeEdge(RemovedEdges.back());
     545            ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
     546          }
    560547        }
    561         Vector PositionUpdate = stepwidth * currentGradient;
    562         LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
    563 
    564         // add PositionUpdate for all nodes in the bondside_set
    565         for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set[side].begin();
    566             setiter != bondside_set[side].end(); ++setiter) {
     548        // perform limited-horizon BFS
     549        BoostGraphHelpers::Nodeset_t bondside_set;
     550        BreadthFirstSearchGatherer::distance_map_t distance_map;
     551        bondside_set = NodeGatherer(walker.getId(), max_distance);
     552        distance_map = NodeGatherer.getDistances();
     553        std::sort(bondside_set.begin(), bondside_set.end());
     554
     555        // re-add edge
     556        for (RemovedEdges_t::const_iterator edgeiter = RemovedEdges.begin();
     557            edgeiter != RemovedEdges.end(); ++edgeiter)
     558          BGcreator.addEdge(edgeiter->first, edgeiter->second);
     559
     560        // update position with dampening factor on the discovered bonds
     561        for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set.begin();
     562            setiter != bondside_set.end(); ++setiter) {
    567563          const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
    568             = distance_map[side].find(*setiter);
    569           ASSERT( diter != distance_map[side].end(),
     564            = distance_map.find(*setiter);
     565          ASSERT( diter != distance_map.end(),
    570566              "ForceAnnealing() - could not find distance to an atom.");
    571567          const double factor = pow(damping_factor, diter->second+1);
     
    581577          }
    582578        }
     579      } else {
     580        // simple atomic annealing, i.e. damping factor of 1
     581        LOG(3, "DEBUG: Update for atom #" << walker.getId() << " will be " << PositionUpdate);
     582        GatheredUpdates.insert(
     583            std::make_pair(
     584                walker.getId(),
     585                PositionUpdate) );
    583586      }
    584587    }
     
    593596    }
    594597
    595     // remove center of weight translation from gathered updates
    596     Vector CommonTranslation;
    597     for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
    598         iter != GatheredUpdates.end(); ++iter) {
    599       const Vector &update = iter->second;
    600       CommonTranslation += update;
    601     }
    602     CommonTranslation *= 1./(double)GatheredUpdates.size();
    603     LOG(3, "DEBUG: Subtracting common translation " << CommonTranslation
    604         << " from all updates.");
     598//    // remove center of weight translation from gathered updates
     599//    Vector CommonTranslation;
     600//    for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
     601//        iter != GatheredUpdates.end(); ++iter) {
     602//      const Vector &update = iter->second;
     603//      CommonTranslation += update;
     604//    }
     605//    CommonTranslation *= 1./(double)GatheredUpdates.size();
     606//    LOG(3, "DEBUG: Subtracting common translation " << CommonTranslation
     607//        << " from all updates.");
    605608
    606609    // apply the gathered updates and set remnant gradients for atomic annealing
     
    615618          << ", namely " << *walker);
    616619      walker->setPositionAtStep(_TimeStep,
    617           walker->getPositionAtStep(CurrentTimeStep)
    618           + update - CommonTranslation);
    619 //      walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
     620          walker->getPositionAtStep(CurrentTimeStep) + update); // - CommonTranslation);
    620621    }
    621622
Note: See TracChangeset for help on using the changeset viewer.