source: src/Dynamics/ForceAnnealing.hpp@ 698308

/*
 * ForceAnnealing.hpp
 *
 *  Created on: Aug 02, 2014
 *      Author: heber
 */

#ifndef FORCEANNEALING_HPP_
#define FORCEANNEALING_HPP_

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <algorithm>
#include <functional>
#include <iterator>

#include <boost/bind.hpp>

#include "Atom/atom.hpp"
#include "Atom/AtomSet.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "Dynamics/AtomicForceManipulator.hpp"
#include "Dynamics/BondVectors.hpp"
#include "Fragmentation/ForceMatrix.hpp"
#include "Graph/BoostGraphCreator.hpp"
#include "Graph/BoostGraphHelpers.hpp"
#include "Graph/BreadthFirstSearchGatherer.hpp"
#include "Helpers/helpers.hpp"
#include "Helpers/defs.hpp"
#include "LinearAlgebra/LinearSystemOfEquations.hpp"
#include "LinearAlgebra/MatrixContent.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/VectorContent.hpp"
#include "Thermostats/ThermoStatContainer.hpp"
#include "Thermostats/Thermostat.hpp"
#include "World.hpp"

/** This class is the essential build block for performing structural optimization.
 *
 * Sadly, we have to use some static instances as so far values cannot be passed
 * between actions. Hence, we need to store the current step and the adaptive-
 * step width (we cannot perform a line search, as we have no control over the
 * calculation of the forces).
 *
 * However, we do use the bond graph, i.e. if a single atom needs to be shifted
 * to the left, then the whole molecule left of it is shifted, too. This is
 * controlled by the \a max_distance parameter.
 */
template <class T>
class ForceAnnealing : public AtomicForceManipulator<T>
{
public:
  /** Constructor of class ForceAnnealing.
   *
   * \note We use a fixed delta t of 1.
   *
   * \param _atoms set of atoms to integrate
   * \param _Deltat time step width in atomic units
   * \param _IsAngstroem whether length units are in angstroem or bohr radii
   * \param _maxSteps number of optimization steps to perform
   * \param _max_distance up to this bond order is bond graph taken into account.
   */
  ForceAnnealing(
      AtomSetMixin<T> &_atoms,
      const double _Deltat,
      bool _IsAngstroem,
      const size_t _maxSteps,
      const int _max_distance,
      const double _damping_factor) :
    AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem),
    maxSteps(_maxSteps),
    max_distance(_max_distance),
    damping_factor(_damping_factor)
  {}

  /** Destructor of class ForceAnnealing.
   *
   */
  ~ForceAnnealing()
  {}

  /** Performs Gradient optimization.
   *
   * We assume that forces have just been calculated.
   *
   *
   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
   * \param offset offset in matrix file to the first force component
   * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
   */
  void operator()(
      const int _CurrentTimeStep,
      const size_t _offset,
      const bool _UseBondgraph)
  {
    // make sum of forces equal zero
    AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(_offset, _CurrentTimeStep);

    // are we in initial step? Then set static entities
    Vector maxComponents(zeroVec);
    if (currentStep == 0) {
      currentDeltat = AtomicForceManipulator<T>::Deltat;
      currentStep = 1;
      LOG(2, "DEBUG: Initial step, setting values, current step is #" << currentStep);

      // always use atomic annealing on first step
      anneal(_CurrentTimeStep, _offset, maxComponents);
    } else {
      ++currentStep;
      LOG(2, "DEBUG: current step is #" << currentStep);

      if (_UseBondgraph)
        annealWithBondGraph(_CurrentTimeStep, _offset, maxComponents);
      else
        anneal(_CurrentTimeStep, _offset, maxComponents);
    }

    LOG(1, "STATUS: Largest remaining force components at step #"
        << currentStep << " are " << maxComponents);

    // are we in final step? Remember to reset static entities
    if (currentStep == maxSteps) {
      LOG(2, "DEBUG: Final step, resetting values");
      reset();
    }
  }

  /** Performs Gradient optimization on the atoms.
   *
   * We assume that forces have just been calculated.
   *
   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
   * \param offset offset in matrix file to the first force component
   * \param maxComponents to be filled with maximum force component over all atoms
   */
  void anneal(
      const int CurrentTimeStep,
      const size_t offset,
      Vector &maxComponents)
  {
    bool deltat_decreased = false;
                for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
                    iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
                        // atom's force vector gives steepest descent direction
                        const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
                        const Vector currentPosition = (*iter)->getPosition();
                        const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
                        const Vector currentGradient = (*iter)->getAtomicForce();
      LOG(4, "DEBUG: oldPosition for atom " << **iter << " is " << oldPosition);
      LOG(4, "DEBUG: currentPosition for atom " << **iter << " is " << currentPosition);
      LOG(4, "DEBUG: oldGradient for atom " << **iter << " is " << oldGradient);
      LOG(4, "DEBUG: currentGradient for atom " << **iter << " is " << currentGradient);
//                      LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);

      // we use Barzilai-Borwein update with position reversed to get descent
      const Vector PositionDifference = currentPosition - oldPosition;
                        const Vector GradientDifference = (currentGradient - oldGradient);
                        double stepwidth = 0.;
                        if (GradientDifference.Norm() > MYEPSILON)
        stepwidth = fabs(PositionDifference.ScalarProduct(GradientDifference))/
            GradientDifference.NormSquared();
                        if (fabs(stepwidth) < 1e-10) {
                          // dont' warn in first step, deltat usage normal
                          if (currentStep != 1)
          ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
                          stepwidth = currentDeltat;
                        }
      Vector PositionUpdate = stepwidth * currentGradient;
                        LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);

                        // extract largest components for showing progress of annealing
                        for(size_t i=0;i<NDIM;++i)
                          if (currentGradient[i] > maxComponents[i])
                            maxComponents[i] = currentGradient[i];

                        // steps may go back and forth again (updates are of same magnitude but
                        // have different sign: Check whether this is the case and one step with
                        // deltat to interrupt this sequence
                        if ((currentStep > 1) && (!PositionDifference.IsZero()))
                                if ((PositionUpdate.ScalarProduct(PositionDifference) < 0)
                                    && (fabs(PositionUpdate.NormSquared()-PositionDifference.NormSquared()) < 1e-3)) {
                                  // for convergence we want a null sequence here, too
                                  if (!deltat_decreased) {
            deltat_decreased = true;
            currentDeltat = .5*currentDeltat;
                                  }
                                  LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate
                                      << " > " << PositionDifference
                                      << ", using deltat: " << currentDeltat);
                                  PositionUpdate = currentDeltat * currentGradient;
                        }

                        // finally set new values
                        (*iter)->setPosition(currentPosition + PositionUpdate);
                }
  }

  /** Performs Gradient optimization on the bonds.
   *
   * We assume that forces have just been calculated. These forces are projected
   * onto the bonds and these are annealed subsequently by moving atoms in the
   * bond neighborhood on either side conjunctively.
   *
   *
   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
   * \param offset offset in matrix file to the first force component
   * \param maxComponents to be filled with maximum force component over all atoms
   */
  void annealWithBondGraph(
      const int CurrentTimeStep,
      const size_t offset,
      Vector &maxComponents)
  {
    // get nodes on either side of selected bond via BFS discovery
    BoostGraphCreator BGcreator;
    BGcreator.createFromRange(
        AtomicForceManipulator<T>::atoms.begin(),
        AtomicForceManipulator<T>::atoms.end(),
        AtomicForceManipulator<T>::atoms.size(),
        BreadthFirstSearchGatherer::AlwaysTruePredicate);
    BreadthFirstSearchGatherer NodeGatherer(BGcreator);

    /// We assume that a force is local, i.e. a bond is too short yet and hence
    /// the atom needs to be moved. However, all the adjacent (bound) atoms might
    /// already be at the perfect distance. If we just move the atom alone, we ruin
    /// all the other bonds. Hence, it would be sensible to move every atom found
    /// through the bond graph in the direction of the force as well by the same
    /// PositionUpdate. This is almost what we are going to do.

    /// One issue is: If we need to shorten bond, then we use the PositionUpdate
    /// also on the the other bond partner already. This is because it is in the
    /// direction of the bond. Therefore, the update is actually performed twice on
    /// each bond partner, i.e. the step size is twice as large as it should be.
    /// This problem only occurs when bonds need to be shortened, not when they
    /// need to be made longer (then the force vector is facing the other
    /// direction than the bond vector).
    /// As a remedy we need to average the force on either end of the bond and
    /// check whether each gradient points inwards out or outwards with respect
    /// to the bond and then shift accordingly.

    /// One more issue is that the projection onto the bond directions does not
    /// recover the gradient but may be larger as the bond directions are a
    /// generating system and not a basis (e.g. 3 bonds on a plane where 2 would
    /// suffice to span the plane). To this end, we need to account for the
    /// overestimation and obtain a weighting for each bond.

    // initialize helper class for bond vectors using bonds from range of atoms
    BondVectors bv;
    bv.setFromAtomRange< T >(
        AtomicForceManipulator<T>::atoms.begin(),
        AtomicForceManipulator<T>::atoms.end(),
        currentStep);
    const BondVectors::container_t &sorted_bonds = bv.getSorted();

    // knowing the number of bonds in total, we can setup the storage for the
    // projected forces
    enum whichatom_t {
      leftside=0,
      rightside=1,
      MAX_sides
    };
    std::vector< // time step
      std::vector< // which bond side
        std::vector<double> > // over all bonds
          > projected_forces(2); // one for leftatoms, one for rightatoms (and for both time steps)
    for (size_t i=0;i<2;++i) {
      projected_forces[i].resize(MAX_sides);
      for (size_t j=0;j<MAX_sides;++j)
        projected_forces[i][j].resize(sorted_bonds.size(), 0.);
    }

    // for each atom we need to gather weights and then project the gradient
    typedef std::deque<double> weights_t;
    typedef std::map<atomId_t, weights_t > weights_per_atom_t;
    std::vector<weights_per_atom_t> weights_per_atom(2);
    for (size_t timestep = 0; timestep <= 1; ++timestep) {
      // TODO: We have this extra step in between because of DoOutput copying
      // change this making the code easier to understand
      const size_t CurrentStep = CurrentTimeStep-2*timestep >= 0 ? CurrentTimeStep - 2*timestep : 0;

      // get all bond vectors for this time step (from the perspective of the
      // bonds taken from the currentStep)
      const BondVectors::mapped_t bondvectors = bv.getBondVectorsAtStep(CurrentStep);

      for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin();
          iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
        const atom &walker = *(*iter);
        const Vector &walkerGradient = walker.getAtomicForceAtStep(CurrentStep);
        LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
            << walker << " is " << walkerGradient);

        const BondList& ListOfBonds = walker.getListOfBonds();
        if (walkerGradient.Norm() > MYEPSILON) {

          // gather subset of BondVectors for the current atom
          std::vector<Vector> BondVectors;
          for(BondList::const_iterator bonditer = ListOfBonds.begin();
              bonditer != ListOfBonds.end(); ++bonditer) {
            const bond::ptr &current_bond = *bonditer;
            const BondVectors::mapped_t::const_iterator bviter =
                bondvectors.find(current_bond);
            ASSERT( bviter != bondvectors.end(),
                "ForceAnnealing() - cannot find current_bond ?");
            ASSERT( bviter != bondvectors.end(),
                "ForceAnnealing - cannot find current bond "+toString(*current_bond)
                +" in bonds.");
            BondVectors.push_back(bviter->second);
          }
          LOG(4, "DEBUG: BondVectors for atom #" << walker.getId() << ": " << BondVectors);

          // go through all its bonds and calculate what magnitude is represented
          // by the others i.e. sum of scalar products against other bonds
          std::pair<weights_per_atom_t::iterator, bool> inserter =
              weights_per_atom[timestep].insert(
                  std::make_pair(walker.getId(), weights_t()) );
          ASSERT( inserter.second,
              "ForceAnnealing::operator() - weight map for atom "+toString(walker)
              +" and time step "+toString(timestep)+" already filled?");
          weights_t &weights = inserter.first->second;
          for (std::vector<Vector>::const_iterator iter = BondVectors.begin();
              iter != BondVectors.end(); ++iter) {
            std::vector<double> scps;
            scps.reserve(BondVectors.size());
            std::transform(
                BondVectors.begin(), BondVectors.end(),
                std::back_inserter(scps),
                boost::bind(static_cast< double (*)(double) >(&fabs),
                    boost::bind(&Vector::ScalarProduct, boost::cref(*iter), _1))
                );
            const double scp_sum = std::accumulate(scps.begin(), scps.end(), 0.);
            ASSERT( (scp_sum-1.) > -MYEPSILON,
                "ForceAnnealing() - sum of weights must be equal or larger one but is "
                +toString(scp_sum));
            weights.push_back( 1./scp_sum );
          }
          LOG(4, "DEBUG: Weights for atom #" << walker.getId() << ": " << weights);

          // for testing we check whether all weighted scalar products now come out as 1.
#ifndef NDEBUG
          for (std::vector<Vector>::const_iterator iter = BondVectors.begin();
              iter != BondVectors.end(); ++iter) {
            std::vector<double> scps;
            scps.reserve(BondVectors.size());
            std::transform(
                BondVectors.begin(), BondVectors.end(),
                weights.begin(),
                std::back_inserter(scps),
                boost::bind(static_cast< double (*)(double) >(&fabs),
                    boost::bind(std::multiplies<double>(),
                        boost::bind(&Vector::ScalarProduct, boost::cref(*iter), _1),
                        _2))
                );
            const double scp_sum = std::accumulate(scps.begin(), scps.end(), 0.);
            ASSERT( fabs(scp_sum - 1.) < MYEPSILON,
                "ForceAnnealing::operator() - for BondVector "+toString(*iter)
                +" we have weighted scalar product of "+toString(scp_sum)+" != 1.");
          }
#endif

          // projected gradient over all bonds and place in one of projected_forces
          // using the obtained weights
          {
            weights_t::const_iterator weightiter = weights.begin();
            std::vector<Vector>::const_iterator vectoriter = BondVectors.begin();
            Vector forcesum_check;
            for(BondList::const_iterator bonditer = ListOfBonds.begin();
                bonditer != ListOfBonds.end(); ++bonditer, ++weightiter, ++vectoriter) {
              const bond::ptr &current_bond = *bonditer;
              const Vector &BondVector = *vectoriter;

              std::vector<double> &forcelist = (&walker == current_bond->leftatom) ?
                  projected_forces[timestep][leftside] : projected_forces[timestep][rightside];
              const size_t index = bv.getIndexForBond(current_bond);
              ASSERT( index != (size_t)-1,
                  "ForceAnnealing() - could not find bond "+toString(*current_bond)
                  +" in bondvectors");
              forcelist[index] = (*weightiter)*walkerGradient.ScalarProduct(BondVector);
              LOG(4, "DEBUG: BondVector " << BondVector << " receives projected force of "
                  << forcelist[index]);
              forcesum_check += forcelist[index] * BondVector;
            }
            ASSERT( weightiter == weights.end(),
                "ForceAnnealing - weightiter is not at end when it should be.");
            ASSERT( vectoriter == BondVectors.end(),
                "ForceAnnealing - vectoriter is not at end when it should be.");
            LOG(3, "DEBUG: sum of projected forces is " << forcesum_check);
          }

        } else {
          LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
              << MYEPSILON << " for atom " << walker);
          // note that projected_forces is initialized to full length and filled
          // with zeros. Hence, nothing to do here
        }
      }
    }

    // step through each bond and shift the atoms
    std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end

    LOG(3, "DEBUG: current step is " << currentStep << ", given time step is " << CurrentTimeStep);
    const BondVectors::mapped_t bondvectors = bv.getBondVectorsAtStep(currentStep);

    for (BondVectors::container_t::const_iterator bondsiter = sorted_bonds.begin();
        bondsiter != sorted_bonds.end(); ++bondsiter) {
      const bond::ptr &current_bond = *bondsiter;
      const size_t index = bv.getIndexForBond(current_bond);
      const atom* bondatom[MAX_sides] = {
          current_bond->leftatom,
          current_bond->rightatom
      };

      // remove the edge
#ifndef NDEBUG
      const bool status =
#endif
          BGcreator.removeEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
      ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");

      // gather nodes for either atom
      BoostGraphHelpers::Nodeset_t bondside_set[MAX_sides];
      BreadthFirstSearchGatherer::distance_map_t distance_map[MAX_sides];
      for (size_t side=leftside;side<MAX_sides;++side) {
        bondside_set[side] = NodeGatherer(bondatom[side]->getId(), max_distance);
        distance_map[side] = NodeGatherer.getDistances();
        std::sort(bondside_set[side].begin(), bondside_set[side].end());
      }

      // re-add edge
      BGcreator.addEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());

      // do for both leftatom and rightatom of bond
      for (size_t side = leftside; side < MAX_sides; ++side) {
        const double &bondforce = projected_forces[0][side][index];
        const double &oldbondforce = projected_forces[1][side][index];
        const double bondforcedifference = fabs(bondforce - oldbondforce);
        LOG(4, "DEBUG: bondforce for " << (side == leftside ? "left" : "right")
            << " side of bond is " << bondforce);
        LOG(4, "DEBUG: oldbondforce for " << (side == leftside ? "left" : "right")
            << " side of bond is " << oldbondforce);
        // if difference or bondforce itself is zero, do nothing
        if ((fabs(bondforce) < MYEPSILON) || (fabs(bondforcedifference) < MYEPSILON))
          continue;

        // get BondVector to bond
        const BondVectors::mapped_t::const_iterator bviter =
            bondvectors.find(current_bond);
        ASSERT( bviter != bondvectors.end(),
            "ForceAnnealing() - cannot find current_bond ?");
        ASSERT( fabs(bviter->second.Norm() -1.) < MYEPSILON,
            "ForceAnnealing() - norm of BondVector is not one");
        const Vector &BondVector = bviter->second;

        // calculate gradient and position differences for stepwidth
        const Vector currentGradient = bondforce * BondVector;
        LOG(4, "DEBUG: current projected gradient for "
            << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient);
        const Vector &oldPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
        const Vector &currentPosition = bondatom[side]->getPositionAtStep(currentStep);
        const Vector PositionDifference = currentPosition - oldPosition;
        LOG(4, "DEBUG: old position is " << oldPosition);
        LOG(4, "DEBUG: current position is " << currentPosition);
        LOG(4, "DEBUG: difference in position is " << PositionDifference);
        LOG(4, "DEBUG: abs. difference in forces is " << bondforcedifference);
        LOG(4, "DEBUG: bondvector is " << BondVector);

        // calculate step width
        double stepwidth =
            fabs(PositionDifference.ScalarProduct(BondVector))/bondforcedifference;
        if (fabs(stepwidth) < 1e-10) {
          // dont' warn in first step, deltat usage normal
          if (currentStep != 1)
            ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
          stepwidth = currentDeltat;
        }
        Vector PositionUpdate = stepwidth * currentGradient;
        LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);

        // add PositionUpdate for all nodes in the bondside_set
        for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set[side].begin();
            setiter != bondside_set[side].end(); ++setiter) {
          const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
            = distance_map[side].find(*setiter);
          ASSERT( diter != distance_map[side].end(),
              "ForceAnnealing() - could not find distance to an atom.");
          const double factor = pow(damping_factor, diter->second+1);
          LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
              << factor << "*" << PositionUpdate);
          if (GatheredUpdates.count((*setiter))) {
            GatheredUpdates[(*setiter)] += factor*PositionUpdate;
          } else {
            GatheredUpdates.insert(
                std::make_pair(
                    (*setiter),
                    factor*PositionUpdate) );
          }
        }
      }
    }

    for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
      atom &walker = *(*iter);
      // extract largest components for showing progress of annealing
      const Vector currentGradient = walker.getAtomicForce();
      for(size_t i=0;i<NDIM;++i)
        if (currentGradient[i] > maxComponents[i])
          maxComponents[i] = currentGradient[i];
    }

    // apply the gathered updates
    for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
        iter != GatheredUpdates.end(); ++iter) {
      const atomId_t &atomid = iter->first;
      const Vector &update = iter->second;
      atom* const walker = World::getInstance().getAtom(AtomById(atomid));
      ASSERT( walker != NULL,
          "ForceAnnealing() - walker with id "+toString(atomid)+" has suddenly disappeared.");
      LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
          << ", namely " << *walker);
      walker->setPosition( walker->getPosition() + update );
    }
  }

  /** Reset function to unset static entities and artificial velocities.
   *
   */
  void reset()
  {
    currentDeltat = 0.;
    currentStep = 0;
  }

private:
  //!> contains the current step in relation to maxsteps
  static size_t currentStep;
  //!> contains the maximum number of steps, determines initial and final step with currentStep
  size_t maxSteps;
  static double currentDeltat;
  //!> minimum deltat for internal while loop (adaptive step width)
  static double MinimumDeltat;
  //!> contains the maximum bond graph distance up to which shifts of a single atom are spread
  const int max_distance;
  //!> the shifted is dampened by this factor with the power of the bond graph distance to the shift causing atom
  const double damping_factor;
};

template <class T>
double ForceAnnealing<T>::currentDeltat = 0.;
template <class T>
size_t ForceAnnealing<T>::currentStep = 0;
template <class T>
double ForceAnnealing<T>::MinimumDeltat = 1e-8;

#endif /* FORCEANNEALING_HPP_ */