1 | /*
|
---|
2 | * ForceAnnealing.hpp
|
---|
3 | *
|
---|
4 | * Created on: Aug 02, 2014
|
---|
5 | * Author: heber
|
---|
6 | */
|
---|
7 |
|
---|
8 | #ifndef FORCEANNEALING_HPP_
|
---|
9 | #define FORCEANNEALING_HPP_
|
---|
10 |
|
---|
11 | // include config.h
|
---|
12 | #ifdef HAVE_CONFIG_H
|
---|
13 | #include <config.h>
|
---|
14 | #endif
|
---|
15 |
|
---|
16 | #include "Atom/atom.hpp"
|
---|
17 | #include "Atom/AtomSet.hpp"
|
---|
18 | #include "CodePatterns/Assert.hpp"
|
---|
19 | #include "CodePatterns/Info.hpp"
|
---|
20 | #include "CodePatterns/Log.hpp"
|
---|
21 | #include "CodePatterns/Verbose.hpp"
|
---|
22 | #include "Descriptors/AtomIdDescriptor.hpp"
|
---|
23 | #include "Dynamics/AtomicForceManipulator.hpp"
|
---|
24 | #include "Fragmentation/ForceMatrix.hpp"
|
---|
25 | #include "Graph/BoostGraphCreator.hpp"
|
---|
26 | #include "Graph/BoostGraphHelpers.hpp"
|
---|
27 | #include "Graph/BreadthFirstSearchGatherer.hpp"
|
---|
28 | #include "Helpers/helpers.hpp"
|
---|
29 | #include "Helpers/defs.hpp"
|
---|
30 | #include "LinearAlgebra/Vector.hpp"
|
---|
31 | #include "Thermostats/ThermoStatContainer.hpp"
|
---|
32 | #include "Thermostats/Thermostat.hpp"
|
---|
33 | #include "World.hpp"
|
---|
34 |
|
---|
35 | /** This class is the essential build block for performing structural optimization.
|
---|
36 | *
|
---|
37 | * Sadly, we have to use some static instances as so far values cannot be passed
|
---|
38 | * between actions. Hence, we need to store the current step and the adaptive-
|
---|
39 | * step width (we cannot perform a line search, as we have no control over the
|
---|
40 | * calculation of the forces).
|
---|
41 | *
|
---|
42 | * However, we do use the bond graph, i.e. if a single atom needs to be shifted
|
---|
43 | * to the left, then the whole molecule left of it is shifted, too. This is
|
---|
44 | * controlled by the \a max_distance parameter.
|
---|
45 | */
|
---|
46 | template <class T>
|
---|
47 | class ForceAnnealing : public AtomicForceManipulator<T>
|
---|
48 | {
|
---|
49 | public:
|
---|
50 | /** Constructor of class ForceAnnealing.
|
---|
51 | *
|
---|
52 | * \note We use a fixed delta t of 1.
|
---|
53 | *
|
---|
54 | * \param _atoms set of atoms to integrate
|
---|
55 | * \param _Deltat time step width in atomic units
|
---|
56 | * \param _IsAngstroem whether length units are in angstroem or bohr radii
|
---|
57 | * \param _maxSteps number of optimization steps to perform
|
---|
58 | * \param _max_distance up to this bond order is bond graph taken into account.
|
---|
59 | */
|
---|
60 | ForceAnnealing(
|
---|
61 | AtomSetMixin<T> &_atoms,
|
---|
62 | const double _Deltat,
|
---|
63 | bool _IsAngstroem,
|
---|
64 | const size_t _maxSteps,
|
---|
65 | const int _max_distance,
|
---|
66 | const double _damping_factor) :
|
---|
67 | AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem),
|
---|
68 | maxSteps(_maxSteps),
|
---|
69 | max_distance(_max_distance),
|
---|
70 | damping_factor(_damping_factor)
|
---|
71 | {}
|
---|
72 |
|
---|
73 | /** Destructor of class ForceAnnealing.
|
---|
74 | *
|
---|
75 | */
|
---|
76 | ~ForceAnnealing()
|
---|
77 | {}
|
---|
78 |
|
---|
79 | /** Performs Gradient optimization.
|
---|
80 | *
|
---|
81 | * We assume that forces have just been calculated.
|
---|
82 | *
|
---|
83 | *
|
---|
84 | * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
|
---|
85 | * \param offset offset in matrix file to the first force component
|
---|
86 | * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
|
---|
87 | */
|
---|
88 | void operator()(
|
---|
89 | const int _CurrentTimeStep,
|
---|
90 | const size_t _offset,
|
---|
91 | const bool _UseBondgraph)
|
---|
92 | {
|
---|
93 | // make sum of forces equal zero
|
---|
94 | AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(_offset, _CurrentTimeStep);
|
---|
95 |
|
---|
96 | // are we in initial step? Then set static entities
|
---|
97 | Vector maxComponents(zeroVec);
|
---|
98 | if (currentStep == 0) {
|
---|
99 | currentDeltat = AtomicForceManipulator<T>::Deltat;
|
---|
100 | currentStep = 1;
|
---|
101 | LOG(2, "DEBUG: Initial step, setting values, current step is #" << currentStep);
|
---|
102 |
|
---|
103 | // always use atomic annealing on first step
|
---|
104 | anneal(_CurrentTimeStep, _offset, maxComponents);
|
---|
105 | } else {
|
---|
106 | ++currentStep;
|
---|
107 | LOG(2, "DEBUG: current step is #" << currentStep);
|
---|
108 |
|
---|
109 | if (_UseBondgraph)
|
---|
110 | annealWithBondGraph(_CurrentTimeStep, _offset, maxComponents);
|
---|
111 | else
|
---|
112 | anneal(_CurrentTimeStep, _offset, maxComponents);
|
---|
113 | }
|
---|
114 |
|
---|
115 | LOG(1, "STATUS: Largest remaining force components at step #"
|
---|
116 | << currentStep << " are " << maxComponents);
|
---|
117 |
|
---|
118 | // are we in final step? Remember to reset static entities
|
---|
119 | if (currentStep == maxSteps) {
|
---|
120 | LOG(2, "DEBUG: Final step, resetting values");
|
---|
121 | reset();
|
---|
122 | }
|
---|
123 | }
|
---|
124 |
|
---|
125 | /** Helper function to calculate the Barzilai-Borwein stepwidth.
|
---|
126 | *
|
---|
127 | * \param _PositionDifference difference in position between current and last step
|
---|
128 | * \param _GradientDifference difference in gradient between current and last step
|
---|
129 | * \return step width according to Barzilai-Borwein
|
---|
130 | */
|
---|
131 | double getBarzilaiBorweinStepwidth(const Vector &_PositionDifference, const Vector &_GradientDifference)
|
---|
132 | {
|
---|
133 | double stepwidth = 0.;
|
---|
134 | if (_GradientDifference.NormSquared() > MYEPSILON)
|
---|
135 | stepwidth = fabs(_PositionDifference.ScalarProduct(_GradientDifference))/
|
---|
136 | _GradientDifference.NormSquared();
|
---|
137 | if (fabs(stepwidth) < 1e-10) {
|
---|
138 | // dont' warn in first step, deltat usage normal
|
---|
139 | if (currentStep != 1)
|
---|
140 | ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
|
---|
141 | stepwidth = currentDeltat;
|
---|
142 | }
|
---|
143 | return stepwidth;
|
---|
144 | }
|
---|
145 |
|
---|
146 | /** Performs Gradient optimization on the atoms.
|
---|
147 | *
|
---|
148 | * We assume that forces have just been calculated.
|
---|
149 | *
|
---|
150 | * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
|
---|
151 | * \param offset offset in matrix file to the first force component
|
---|
152 | * \param maxComponents to be filled with maximum force component over all atoms
|
---|
153 | */
|
---|
154 | void anneal(
|
---|
155 | const int CurrentTimeStep,
|
---|
156 | const size_t offset,
|
---|
157 | Vector &maxComponents)
|
---|
158 | {
|
---|
159 | for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
|
---|
160 | iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
|
---|
161 | // atom's force vector gives steepest descent direction
|
---|
162 | const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
|
---|
163 | const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
|
---|
164 | const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
|
---|
165 | const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
|
---|
166 | LOG(4, "DEBUG: oldPosition for atom " << **iter << " is " << oldPosition);
|
---|
167 | LOG(4, "DEBUG: currentPosition for atom " << **iter << " is " << currentPosition);
|
---|
168 | LOG(4, "DEBUG: oldGradient for atom " << **iter << " is " << oldGradient);
|
---|
169 | LOG(4, "DEBUG: currentGradient for atom " << **iter << " is " << currentGradient);
|
---|
170 | // LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
|
---|
171 |
|
---|
172 | // we use Barzilai-Borwein update with position reversed to get descent
|
---|
173 | const double stepwidth = getBarzilaiBorweinStepwidth(
|
---|
174 | currentPosition - oldPosition, currentGradient - oldGradient);
|
---|
175 | Vector PositionUpdate = stepwidth * currentGradient;
|
---|
176 | LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
|
---|
177 |
|
---|
178 | // extract largest components for showing progress of annealing
|
---|
179 | for(size_t i=0;i<NDIM;++i)
|
---|
180 | maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
|
---|
181 |
|
---|
182 | // are we in initial step? Then don't check against velocity
|
---|
183 | if ((currentStep > 1) && (!(*iter)->getAtomicVelocity().IsZero()))
|
---|
184 | // update with currentDelta tells us how the current gradient relates to
|
---|
185 | // the last one: If it has become larger, reduce currentDelta
|
---|
186 | if ((PositionUpdate.ScalarProduct((*iter)->getAtomicVelocity()) < 0)
|
---|
187 | && (currentDeltat > MinimumDeltat)) {
|
---|
188 | currentDeltat = .5*currentDeltat;
|
---|
189 | LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate.NormSquared()
|
---|
190 | << " > " << (*iter)->getAtomicVelocity().NormSquared()
|
---|
191 | << ", decreasing deltat: " << currentDeltat);
|
---|
192 | PositionUpdate = currentDeltat * currentGradient;
|
---|
193 | }
|
---|
194 | // finally set new values
|
---|
195 | (*iter)->setPosition(currentPosition + PositionUpdate);
|
---|
196 | (*iter)->setAtomicVelocity(PositionUpdate);
|
---|
197 | //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
|
---|
198 | // (*iter)->VelocityVerletUpdateU((*iter)->getId(), CurrentTimeStep-1, Deltat, IsAngstroem);
|
---|
199 | }
|
---|
200 | }
|
---|
201 |
|
---|
202 | /** Performs Gradient optimization on the bonds.
|
---|
203 | *
|
---|
204 | * We assume that forces have just been calculated. These forces are projected
|
---|
205 | * onto the bonds and these are annealed subsequently by moving atoms in the
|
---|
206 | * bond neighborhood on either side conjunctively.
|
---|
207 | *
|
---|
208 | *
|
---|
209 | * \param CurrentTimeStep current time step (i.e. t where \f$ t + \Delta t \f$ is in the sense of the velocity verlet)
|
---|
210 | * \param offset offset in matrix file to the first force component
|
---|
211 | * \param maxComponents to be filled with maximum force component over all atoms
|
---|
212 | */
|
---|
213 | void annealWithBondGraph(
|
---|
214 | const int CurrentTimeStep,
|
---|
215 | const size_t offset,
|
---|
216 | Vector &maxComponents)
|
---|
217 | {
|
---|
218 | // get nodes on either side of selected bond via BFS discovery
|
---|
219 | // std::vector<atomId_t> atomids;
|
---|
220 | // for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
|
---|
221 | // iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
|
---|
222 | // atomids.push_back((*iter)->getId());
|
---|
223 | // }
|
---|
224 | // ASSERT( atomids.size() == AtomicForceManipulator<T>::atoms.size(),
|
---|
225 | // "ForceAnnealing() - could not gather all atomic ids?");
|
---|
226 | BoostGraphCreator BGcreator;
|
---|
227 | BGcreator.createFromRange(
|
---|
228 | AtomicForceManipulator<T>::atoms.begin(),
|
---|
229 | AtomicForceManipulator<T>::atoms.end(),
|
---|
230 | AtomicForceManipulator<T>::atoms.size(),
|
---|
231 | BreadthFirstSearchGatherer::AlwaysTruePredicate);
|
---|
232 | BreadthFirstSearchGatherer NodeGatherer(BGcreator);
|
---|
233 |
|
---|
234 | std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
|
---|
235 | for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
|
---|
236 | iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
|
---|
237 | // atom's force vector gives steepest descent direction
|
---|
238 | const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
|
---|
239 | const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
|
---|
240 | const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
|
---|
241 | const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
|
---|
242 | LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
|
---|
243 |
|
---|
244 | // we use Barzilai-Borwein update with position reversed to get descent
|
---|
245 | const double stepwidth = getBarzilaiBorweinStepwidth(
|
---|
246 | currentPosition - oldPosition, currentGradient - oldGradient);
|
---|
247 | const Vector PositionUpdate = stepwidth * currentGradient;
|
---|
248 | LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
|
---|
249 |
|
---|
250 | // // add update to central atom
|
---|
251 | // const atomId_t atomid = (*iter)->getId();
|
---|
252 | // if (GatheredUpdates.count(atomid)) {
|
---|
253 | // GatheredUpdates[atomid] += PositionUpdate;
|
---|
254 | // } else
|
---|
255 | // GatheredUpdates.insert( std::make_pair(atomid, PositionUpdate) );
|
---|
256 |
|
---|
257 | // We assume that a force is local, i.e. a bond is too short yet and hence
|
---|
258 | // the atom needs to be moved. However, all the adjacent (bound) atoms might
|
---|
259 | // already be at the perfect distance. If we just move the atom alone, we ruin
|
---|
260 | // all the other bonds. Hence, it would be sensible to move every atom found
|
---|
261 | // through the bond graph in the direction of the force as well by the same
|
---|
262 | // PositionUpdate. This is just what we are going to do.
|
---|
263 |
|
---|
264 | /// get all nodes from bonds in the direction of the current force
|
---|
265 |
|
---|
266 | // remove edges facing in the wrong direction
|
---|
267 | std::vector<bond::ptr> removed_bonds;
|
---|
268 | const BondList& ListOfBonds = (*iter)->getListOfBonds();
|
---|
269 | for(BondList::const_iterator bonditer = ListOfBonds.begin();
|
---|
270 | bonditer != ListOfBonds.end(); ++bonditer) {
|
---|
271 | const bond ¤t_bond = *(*bonditer);
|
---|
272 | LOG(2, "DEBUG: Looking at bond " << current_bond);
|
---|
273 | Vector BondVector = (*iter)->getPositionAtStep(CurrentTimeStep);
|
---|
274 | BondVector -= ((*iter)->getId() == current_bond.rightatom->getId())
|
---|
275 | ? current_bond.rightatom->getPositionAtStep(CurrentTimeStep) : current_bond.leftatom->getPositionAtStep(CurrentTimeStep);
|
---|
276 | BondVector.Normalize();
|
---|
277 | if (BondVector.ScalarProduct(currentGradient) < 0) {
|
---|
278 | removed_bonds.push_back(*bonditer);
|
---|
279 | #ifndef NDEBUG
|
---|
280 | const bool status =
|
---|
281 | #endif
|
---|
282 | BGcreator.removeEdge(current_bond.leftatom->getId(), current_bond.rightatom->getId());
|
---|
283 | ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
|
---|
284 | }
|
---|
285 | }
|
---|
286 | BoostGraphHelpers::Nodeset_t bondside_set = NodeGatherer((*iter)->getId(), max_distance);
|
---|
287 | const BreadthFirstSearchGatherer::distance_map_t &distance_map = NodeGatherer.getDistances();
|
---|
288 | std::sort(bondside_set.begin(), bondside_set.end());
|
---|
289 | // re-add those edges
|
---|
290 | for (std::vector<bond::ptr>::const_iterator bonditer = removed_bonds.begin();
|
---|
291 | bonditer != removed_bonds.end(); ++bonditer)
|
---|
292 | BGcreator.addEdge((*bonditer)->leftatom->getId(), (*bonditer)->rightatom->getId());
|
---|
293 |
|
---|
294 | // apply PositionUpdate to all nodes in the bondside_set
|
---|
295 | for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set.begin();
|
---|
296 | setiter != bondside_set.end(); ++setiter) {
|
---|
297 | const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
|
---|
298 | = distance_map.find(*setiter);
|
---|
299 | ASSERT( diter != distance_map.end(),
|
---|
300 | "ForceAnnealing() - could not find distance to an atom.");
|
---|
301 | const double factor = pow(damping_factor, diter->second);
|
---|
302 | LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
|
---|
303 | << factor << "*" << PositionUpdate);
|
---|
304 | if (GatheredUpdates.count((*setiter))) {
|
---|
305 | GatheredUpdates[(*setiter)] += factor*PositionUpdate;
|
---|
306 | } else {
|
---|
307 | GatheredUpdates.insert(
|
---|
308 | std::make_pair(
|
---|
309 | (*setiter),
|
---|
310 | factor*PositionUpdate) );
|
---|
311 | }
|
---|
312 | }
|
---|
313 |
|
---|
314 | // extract largest components for showing progress of annealing
|
---|
315 | for(size_t i=0;i<NDIM;++i)
|
---|
316 | maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
|
---|
317 | }
|
---|
318 | // apply the gathered updates
|
---|
319 | for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
|
---|
320 | iter != GatheredUpdates.end(); ++iter) {
|
---|
321 | const atomId_t &atomid = iter->first;
|
---|
322 | const Vector &update = iter->second;
|
---|
323 | atom* const walker = World::getInstance().getAtom(AtomById(atomid));
|
---|
324 | ASSERT( walker != NULL,
|
---|
325 | "ForceAnnealing() - walker with id "+toString(atomid)+" has suddenly disappeared.");
|
---|
326 | LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
|
---|
327 | << ", namely " << *walker);
|
---|
328 | walker->setPosition( walker->getPosition() + update );
|
---|
329 | }
|
---|
330 | }
|
---|
331 |
|
---|
332 | /** Reset function to unset static entities and artificial velocities.
|
---|
333 | *
|
---|
334 | */
|
---|
335 | void reset()
|
---|
336 | {
|
---|
337 | currentDeltat = 0.;
|
---|
338 | currentStep = 0;
|
---|
339 | }
|
---|
340 |
|
---|
341 | private:
|
---|
342 | //!> contains the current step in relation to maxsteps
|
---|
343 | static size_t currentStep;
|
---|
344 | //!> contains the maximum number of steps, determines initial and final step with currentStep
|
---|
345 | size_t maxSteps;
|
---|
346 | static double currentDeltat;
|
---|
347 | //!> minimum deltat for internal while loop (adaptive step width)
|
---|
348 | static double MinimumDeltat;
|
---|
349 | //!> contains the maximum bond graph distance up to which shifts of a single atom are spread
|
---|
350 | const int max_distance;
|
---|
351 | //!> the shifted is dampened by this factor with the power of the bond graph distance to the shift causing atom
|
---|
352 | const double damping_factor;
|
---|
353 | };
|
---|
354 |
|
---|
355 | template <class T>
|
---|
356 | double ForceAnnealing<T>::currentDeltat = 0.;
|
---|
357 | template <class T>
|
---|
358 | size_t ForceAnnealing<T>::currentStep = 0;
|
---|
359 | template <class T>
|
---|
360 | double ForceAnnealing<T>::MinimumDeltat = 1e-8;
|
---|
361 |
|
---|
362 | #endif /* FORCEANNEALING_HPP_ */
|
---|