1 | /*
|
---|
2 | * MinimiseConstrainedPotential.hpp
|
---|
3 | *
|
---|
4 | * Created on: Feb 23, 2011
|
---|
5 | * Author: heber
|
---|
6 | */
|
---|
7 |
|
---|
8 | #ifndef MINIMISECONSTRAINEDPOTENTIAL_HPP_
|
---|
9 | #define MINIMISECONSTRAINEDPOTENTIAL_HPP_
|
---|
10 |
|
---|
11 | // include config.h
|
---|
12 | #ifdef HAVE_CONFIG_H
|
---|
13 | #include <config.h>
|
---|
14 | #endif
|
---|
15 |
|
---|
16 | class atom;
|
---|
17 |
|
---|
18 | #include <vector>
|
---|
19 | #include <map>
|
---|
20 |
|
---|
21 | #include "World.hpp"
|
---|
22 |
|
---|
23 | /** Structure to contain parameters needed for evaluation of constraint potential.
|
---|
24 | *
|
---|
25 | */
|
---|
26 | class MinimiseConstrainedPotential
|
---|
27 | {
|
---|
28 | public:
|
---|
29 | /** Constructor.
|
---|
30 | *
|
---|
31 | * @param _atoms set of atoms to operate on
|
---|
32 | * \param _PermutationMap on return: mapping between the atom label of the initial and the final configuration
|
---|
33 | * @return
|
---|
34 | */
|
---|
35 | MinimiseConstrainedPotential(World::AtomComposite &_atoms, std::map<atom*, atom *> &_PermutationMap);
|
---|
36 |
|
---|
37 | /** Destructor.
|
---|
38 | *
|
---|
39 | * @return
|
---|
40 | */
|
---|
41 | ~MinimiseConstrainedPotential();
|
---|
42 |
|
---|
43 | /** Minimizes the extra potential for constrained molecular dynamics and gives forces and the constrained potential energy.
|
---|
44 | * We do the following:
|
---|
45 | * -# Generate a distance list from all source to all target points
|
---|
46 | * -# Sort this per source point
|
---|
47 | * -# Take for each source point the target point with minimum distance, use this as initial permutation
|
---|
48 | * -# check whether molecule::ConstrainedPotential() is greater than injective penalty
|
---|
49 | * -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential.
|
---|
50 | * -# Next, we only apply transformations that keep the injectivity of the permutations list.
|
---|
51 | * -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target
|
---|
52 | * point and try to change it for one with lesser distance, or for the next one with greater distance, but only
|
---|
53 | * if this decreases the conditional potential.
|
---|
54 | * -# finished.
|
---|
55 | * -# Then, we calculate the forces by taking the spatial derivative, where we scale the potential to such a degree,
|
---|
56 | * that the total force is always pointing in direction of the constraint force (ensuring that we move in the
|
---|
57 | * right direction).
|
---|
58 | * -# Finally, we calculate the potential energy and return.
|
---|
59 | * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
|
---|
60 | * \param endstep step giving final position in constrained MD
|
---|
61 | * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
|
---|
62 | * \sa molecule::VerletForceIntegration()
|
---|
63 | * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2)
|
---|
64 | * \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
|
---|
65 | * to ensure they're properties (e.g. constants[2] always greater than the energy of the system).
|
---|
66 | * \bug this all is not O(N log N) but O(N^2)
|
---|
67 | */
|
---|
68 | double operator()(int startstep, int endstep, bool IsAngstroem);
|
---|
69 |
|
---|
70 | /** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
|
---|
71 | * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
|
---|
72 | * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
|
---|
73 | */
|
---|
74 | void EvaluateConstrainedForces(ForceMatrix *Force);
|
---|
75 |
|
---|
76 | private:
|
---|
77 | typedef std::pair < double, atom* > DistancePair;
|
---|
78 | typedef std::multimap < double, atom* > DistanceMap;
|
---|
79 | typedef std::pair < DistanceMap::iterator, bool> DistanceTestPair;
|
---|
80 |
|
---|
81 | World::AtomComposite atoms;
|
---|
82 | int startstep; //!< start configuration (MDStep in atom::trajectory)
|
---|
83 | int endstep; //!< end configuration (MDStep in atom::trajectory)
|
---|
84 | 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$ )
|
---|
85 | std::map<atom *, DistanceMap> DistanceList; //!< distance list of each atom to each atom
|
---|
86 | std::map<atom *, DistanceMap::iterator> StepList; //!< iterator to ascend through NearestNeighbours \a **DistanceList
|
---|
87 | std::map<atom *, unsigned int> DoubleList; //!< count of which sources want to move to this target, basically the injective measure (>1 -> not injective)
|
---|
88 | std::map<atom *, DistanceMap::iterator> DistanceIterators; //!< marks which was the last picked target as injective candidate with smallest distance
|
---|
89 | bool IsAngstroem; //!< whether coordinates are in angstroem (true) or bohrradius (false)
|
---|
90 | double *PenaltyConstants; //!< penalty constant in front of each term
|
---|
91 |
|
---|
92 | /** \f$O(N^2)\f$ operation of calculation distance between each atom pair and putting into DistanceList.
|
---|
93 | */
|
---|
94 | void FillDistanceList();
|
---|
95 |
|
---|
96 | /** Initialize lists.
|
---|
97 | */
|
---|
98 | void CreateInitialLists();
|
---|
99 |
|
---|
100 | /** Permutes \a **&PermutationMap until the penalty is below constants[2].
|
---|
101 | */
|
---|
102 | void MakeInjectivePermutation();
|
---|
103 |
|
---|
104 | /** Calculates the number of doubles in PermutationMap.
|
---|
105 | */
|
---|
106 | unsigned int CalculateDoubleList();
|
---|
107 |
|
---|
108 | /** Print the current permutation map.
|
---|
109 | */
|
---|
110 | void PrintPermutationMap() const;
|
---|
111 |
|
---|
112 | /** Evaluates the potential energy used for constrained molecular dynamics.
|
---|
113 | * \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$
|
---|
114 | * 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
|
---|
115 | * trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
|
---|
116 | * Note that for the second term we have to solve the following linear system:
|
---|
117 | * \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,
|
---|
118 | * offset vector \f$p_1\f$ in direction \f$n_1\f$, offset vector \f$p_2\f$ in direction \f$n_2\f$,
|
---|
119 | * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines.
|
---|
120 | * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
|
---|
121 | * \return potential energy
|
---|
122 | * \note This routine is scaling quadratically which is not optimal.
|
---|
123 | * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
|
---|
124 | */
|
---|
125 | double ConstrainedPotential();
|
---|
126 |
|
---|
127 | /** Try the next nearest neighbour in order to make the permutation map injective.
|
---|
128 | * \param *Walker atom to change its target
|
---|
129 | * \param &OldPotential old value of constraint potential to see if we do better with new target
|
---|
130 | */
|
---|
131 | double TryNextNearestNeighbourForInjectivePermutation(atom *Walker, double &OldPotential);
|
---|
132 |
|
---|
133 | /** Penalizes atoms heading to same target.
|
---|
134 | * \param *Walker atom to check against others
|
---|
135 | * \return \a penalty times the number of equal targets
|
---|
136 | */
|
---|
137 | double PenalizeEqualTargets(atom *Walker);
|
---|
138 |
|
---|
139 | /** Penalizes long trajectories.
|
---|
140 | * \param *Walker atom to check against others
|
---|
141 | * \return penalty times each distance
|
---|
142 | */
|
---|
143 | double SumDistanceOfTrajectories(atom *Walker);
|
---|
144 |
|
---|
145 | };
|
---|
146 |
|
---|
147 |
|
---|
148 | #endif /* MINIMISECONSTRAINEDPOTENTIAL_HPP_ */
|
---|