source: src/Dynamics/MinimiseConstrainedPotential.hpp@ 47d041

/*
 * MinimiseConstrainedPotential.hpp
 *
 *  Created on: Feb 23, 2011
 *      Author: heber
 */

#ifndef MINIMISECONSTRAINEDPOTENTIAL_HPP_
#define MINIMISECONSTRAINEDPOTENTIAL_HPP_

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

class atom;

#include <vector>
#include <map>

#include "molecule.hpp"

/** Structure to contain parameters needed for evaluation of constraint potential.
 *
 */
class MinimiseConstrainedPotential
{
public:
  /** Constructor.
   *
   * @param _atoms set of atoms to operate on
   * \param _PermutationMap on return: mapping between the atom label of the initial and the final configuration
   * @return
   */
  MinimiseConstrainedPotential(molecule::atomSet &_atoms, std::map<atom*, atom *> &_PermutationMap);

  /** Destructor.
   *
   * @return
   */
  ~MinimiseConstrainedPotential();

  /** Minimises the extra potential for constrained molecular dynamics and gives forces and the constrained potential energy.
   * We do the following:
   *  -# Generate a distance list from all source to all target points
   *  -# Sort this per source point
   *  -# Take for each source point the target point with minimum distance, use this as initial permutation
   *  -# check whether molecule::ConstrainedPotential() is greater than injective penalty
   *     -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential.
   *  -# Next, we only apply transformations that keep the injectivity of the permutations list.
   *  -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target
   *     point and try to change it for one with lesser distance, or for the next one with greater distance, but only
   *     if this decreases the conditional potential.
   *  -# finished.
   *  -# Then, we calculate the forces by taking the spatial derivative, where we scale the potential to such a degree,
   *     that the total force is always pointing in direction of the constraint force (ensuring that we move in the
   *     right direction).
   *  -# Finally, we calculate the potential energy and return.
   * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
   * \param endstep step giving final position in constrained MD
   * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
   * \sa molecule::VerletForceIntegration()
   * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2)
   * \todo The constrained potential's constants are set to fixed values right now, but they should scale based on checks of the system in order
   *       to ensure they're properties (e.g. constants[2] always greater than the energy of the system).
   * \bug this all is not O(N log N) but O(N^2)
   */
  double operator()(int startstep, int endstep, bool IsAngstroem);

  /** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
   * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
   * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
   */
  void EvaluateConstrainedForces(ForceMatrix *Force);

private:
  typedef std::pair < double, atom* > DistancePair;
  typedef std::multimap < double, atom* > DistanceMap;
  typedef std::pair < DistanceMap::iterator, bool> DistanceTestPair;

  molecule::atomSet atoms;
  int startstep; //!< start configuration (MDStep in atom::trajectory)
  int endstep; //!< end configuration (MDStep in atom::trajectory)
  std::map<atom*, atom *> &PermutationMap; //!< gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
  std::map<atom *, DistanceMap> DistanceList; //!< distance list of each atom to each atom
  std::map<atom *, DistanceMap::iterator> StepList; //!< iterator to ascend through NearestNeighbours \a **DistanceList
  std::map<atom *, unsigned int> DoubleList; //!< count of which sources want to move to this target, basically the injective measure (>1 -> not injective)
  std::map<atom *, DistanceMap::iterator> DistanceIterators; //!< marks which was the last picked target as injective candidate with smallest distance
  bool IsAngstroem; //!< whether coordinates are in angstroem (true) or bohrradius (false)
  double *PenaltyConstants; //!<  penalty constant in front of each term

  /** \f$O(N^2)\f$ operation of calculation distance between each atom pair and putting into DistanceList.
   */
  void FillDistanceList();

  /** Initialize lists.
   */
  void CreateInitialLists();

  /** Permutes \a **&PermutationMap until the penalty is below constants[2].
   */
  void MakeInjectivePermutation();

  /** Calculates the number of doubles in PermutationMap.
   */
  unsigned int CalculateDoubleList();

  /** Print the current permutation map.
   */
  void PrintPermutationMap() const;

  /** Evaluates the potential energy used for constrained molecular dynamics.
   * \f$V_i^{con} = c^{bond} \cdot | r_{P(i)} - R_i | + sum_{i \neq j} C^{min} \cdot \frac{1}{C_{ij}} + C^{inj} \Bigl (1 - \theta \bigl (\prod_{i \neq j} (P(i) - P(j)) \bigr ) \Bigr )\f$
   *     where the first term points to the target in minimum distance, the second is a penalty for trajectories lying too close to each other (\f$C_{ij}\f$ is minimum distance between
   *     trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
   * Note that for the second term we have to solve the following linear system:
   * \f$-c_1 \cdot n_1 + c_2 \cdot n_2 + C \cdot n_3 = - p_2 + p_1\f$, where \f$c_1\f$, \f$c_2\f$ and \f$C\f$ are constants,
   * offset vector \f$p_1\f$ in direction \f$n_1\f$, offset vector \f$p_2\f$ in direction \f$n_2\f$,
   * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines.
   * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
   * \return potential energy
   * \note This routine is scaling quadratically which is not optimal.
   * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
   */
  double ConstrainedPotential();

  /** Try the next nearest neighbour in order to make the permutation map injective.
   * \param *Walker atom to change its target
   * \param &OldPotential old value of constraint potential to see if we do better with new target
   */
  double TryNextNearestNeighbourForInjectivePermutation(atom *Walker, double &OldPotential);

  /** Penalizes atoms heading to same target.
   * \param *Walker atom to check against others
   * \return \a penalty times the number of equal targets
   */
  double PenalizeEqualTargets(atom *Walker);

  /** Penalizes long trajectories.
   * \param *Walker atom to check against others
   * \return penalty times each distance
   */
  double SumDistanceOfTrajectories(atom *Walker);

};


#endif /* MINIMISECONSTRAINEDPOTENTIAL_HPP_ */