source: src/Dynamics/ForceAnnealing.hpp@ bd19c1

Candidate_v1.6.1 ChemicalSpaceEvaluator Exclude_Hydrogens_annealWithBondGraph ForceAnnealing_with_BondGraph_contraction-expansion
Last change on this file since bd19c1 was bd19c1, checked in by Frederik Heber <frederik.heber@…>, 7 years ago

Capping BB update step at 0.2 and general BB stepwidth at maximum of 1.

  • TESTFIX: Python regression tests on ForceAnnealing for five carbon have been modified. With bondgraph results are not as good as before.
  • Property mode set to 100644
File size: 33.3 KB
Line 
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 <algorithm>
17#include <functional>
18#include <iterator>
19#include <math.h>
20
21#include <boost/bind.hpp>
22
23#include "Atom/atom.hpp"
24#include "Atom/AtomSet.hpp"
25#include "CodePatterns/Assert.hpp"
26#include "CodePatterns/Info.hpp"
27#include "CodePatterns/Log.hpp"
28#include "CodePatterns/Verbose.hpp"
29#include "Descriptors/AtomIdDescriptor.hpp"
30#include "Dynamics/AtomicForceManipulator.hpp"
31#include "Dynamics/BondVectors.hpp"
32#include "Fragmentation/ForceMatrix.hpp"
33#include "Graph/BoostGraphCreator.hpp"
34#include "Graph/BoostGraphHelpers.hpp"
35#include "Graph/BreadthFirstSearchGatherer.hpp"
36#include "Helpers/helpers.hpp"
37#include "Helpers/defs.hpp"
38#include "LinearAlgebra/LinearSystemOfEquations.hpp"
39#include "LinearAlgebra/MatrixContent.hpp"
40#include "LinearAlgebra/Vector.hpp"
41#include "LinearAlgebra/VectorContent.hpp"
42#include "Thermostats/ThermoStatContainer.hpp"
43#include "Thermostats/Thermostat.hpp"
44#include "World.hpp"
45
46/** This class is the essential build block for performing structural optimization.
47 *
48 * Sadly, we have to use some static instances as so far values cannot be passed
49 * between actions. Hence, we need to store the current step and the adaptive-
50 * step width (we cannot perform a line search, as we have no control over the
51 * calculation of the forces).
52 *
53 * However, we do use the bond graph, i.e. if a single atom needs to be shifted
54 * to the left, then the whole molecule left of it is shifted, too. This is
55 * controlled by the \a max_distance parameter.
56 */
57template <class T>
58class ForceAnnealing : public AtomicForceManipulator<T>
59{
60public:
61 /** Constructor of class ForceAnnealing.
62 *
63 * \note We use a fixed delta t of 1.
64 *
65 * \param _atoms set of atoms to integrate
66 * \param _Deltat time step width in atomic units
67 * \param _IsAngstroem whether length units are in angstroem or bohr radii
68 * \param _maxSteps number of optimization steps to perform
69 * \param _max_distance up to this bond order is bond graph taken into account.
70 */
71 ForceAnnealing(
72 AtomSetMixin<T> &_atoms,
73 const double _Deltat,
74 bool _IsAngstroem,
75 const size_t _maxSteps,
76 const int _max_distance,
77 const double _damping_factor) :
78 AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem),
79 maxSteps(_maxSteps),
80 max_distance(_max_distance),
81 damping_factor(_damping_factor),
82 FORCE_THRESHOLD(1e-8)
83 {}
84
85 /** Destructor of class ForceAnnealing.
86 *
87 */
88 ~ForceAnnealing()
89 {}
90
91 /** Performs Gradient optimization.
92 *
93 * We assume that forces have just been calculated.
94 *
95 *
96 * \param _TimeStep time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
97 * \param offset offset in matrix file to the first force component
98 * \return false - need to continue annealing, true - may stop because forces very small
99 * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
100 */
101 bool operator()(
102 const int _TimeStep,
103 const size_t _offset,
104 const bool _UseBondgraph)
105 {
106 const int CurrentTimeStep = _TimeStep-1;
107 ASSERT( CurrentTimeStep >= 0,
108 "ForceAnnealing::operator() - a new time step (upon which we work) must already have been copied.");
109
110 // make sum of forces equal zero
111 AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(
112 _offset,
113 CurrentTimeStep);
114
115 // are we in initial step? Then set static entities
116 Vector maxComponents(zeroVec);
117 if (currentStep == 0) {
118 currentDeltat = AtomicForceManipulator<T>::Deltat;
119 currentStep = 1;
120 LOG(2, "DEBUG: Initial step, setting values, current step is #" << currentStep);
121
122 // always use atomic annealing on first step
123 maxComponents = anneal(_TimeStep);
124 } else {
125 ++currentStep;
126 LOG(2, "DEBUG: current step is #" << currentStep);
127
128 // bond graph annealing is always followed by a normal annealing
129 if (_UseBondgraph)
130 maxComponents = annealWithBondGraph_BarzilaiBorwein(_TimeStep);
131 // cannot store RemnantGradient in Atom's Force as it ruins BB stepwidth calculation
132 else
133 maxComponents = anneal_BarzilaiBorwein(_TimeStep);
134 }
135
136
137 LOG(1, "STATUS: Largest remaining force components at step #"
138 << currentStep << " are " << maxComponents);
139
140 // check whether are smaller than threshold
141 bool AnnealingFinished = false;
142 double maxcomp = 0.;
143 for (size_t i=0;i<NDIM;++i)
144 maxcomp = std::max(maxcomp, fabs(maxComponents[i]));
145 if (maxcomp < FORCE_THRESHOLD) {
146 LOG(1, "STATUS: Force components are all less than " << FORCE_THRESHOLD
147 << ", stopping.");
148 currentStep = maxSteps;
149 AnnealingFinished = true;
150 }
151
152 // are we in final step? Remember to reset static entities
153 if (currentStep == maxSteps) {
154 LOG(2, "DEBUG: Final step, resetting values");
155 reset();
156 }
157
158 return AnnealingFinished;
159 }
160
161 /** Helper function to calculate the Barzilai-Borwein stepwidth.
162 *
163 * \param _PositionDifference difference in position between current and last step
164 * \param _GradientDifference difference in gradient between current and last step
165 * \return step width according to Barzilai-Borwein
166 */
167 double getBarzilaiBorweinStepwidth(const Vector &_PositionDifference, const Vector &_GradientDifference)
168 {
169 double stepwidth = 0.;
170 if (_GradientDifference.Norm() > MYEPSILON)
171 stepwidth = fabs(_PositionDifference.ScalarProduct(_GradientDifference))/
172 _GradientDifference.NormSquared();
173 if (fabs(stepwidth) < 1e-10) {
174 // dont' warn in first step, deltat usage normal
175 if (currentStep != 1)
176 ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
177 stepwidth = currentDeltat;
178 }
179 return std::min(1., stepwidth);
180 }
181
182 /** Performs Gradient optimization on the atoms.
183 *
184 * We assume that forces have just been calculated.
185 *
186 * \param _TimeStep time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
187 * \return to be filled with maximum force component over all atoms
188 */
189 Vector anneal(
190 const int _TimeStep)
191 {
192 const int CurrentTimeStep = _TimeStep-1;
193 ASSERT( CurrentTimeStep >= 0,
194 "ForceAnnealing::anneal() - a new time step (upon which we work) must already have been copied.");
195
196 LOG(1, "STATUS: performing simple anneal with default stepwidth " << currentDeltat << " at step #" << currentStep);
197
198 Vector maxComponents;
199 bool deltat_decreased = false;
200 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
201 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
202 // atom's force vector gives steepest descent direction
203 const Vector &currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
204 const Vector &currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
205 LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition);
206 LOG(4, "DEBUG: currentGradient for atom #" << (*iter)->getId() << " is " << currentGradient);
207// LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
208
209 // we use Barzilai-Borwein update with position reversed to get descent
210 double stepwidth = currentDeltat;
211 Vector PositionUpdate = stepwidth * currentGradient;
212 LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
213
214 // extract largest components for showing progress of annealing
215 for(size_t i=0;i<NDIM;++i)
216 maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
217
218 // steps may go back and forth again (updates are of same magnitude but
219 // have different sign: Check whether this is the case and one step with
220 // deltat to interrupt this sequence
221 if (currentStep > 1) {
222 const int OldTimeStep = CurrentTimeStep-1;
223 ASSERT( OldTimeStep >= 0,
224 "ForceAnnealing::anneal() - if currentStep is "+toString(currentStep)
225 +", then there should be at least three time steps.");
226 const Vector &oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
227 const Vector PositionDifference = currentPosition - oldPosition;
228 LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
229 LOG(4, "DEBUG: PositionDifference for atom #" << (*iter)->getId() << " is " << PositionDifference);
230 if ((PositionUpdate.ScalarProduct(PositionDifference) < 0)
231 && (fabs(PositionUpdate.NormSquared()-PositionDifference.NormSquared()) < 1e-3)) {
232 // for convergence we want a null sequence here, too
233 if (!deltat_decreased) {
234 deltat_decreased = true;
235 currentDeltat = .5*currentDeltat;
236 }
237 LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate
238 << " > " << PositionDifference
239 << ", using deltat: " << currentDeltat);
240 PositionUpdate = currentDeltat * currentGradient;
241 }
242 }
243
244 // finally set new values
245 (*iter)->setPositionAtStep(_TimeStep, currentPosition + PositionUpdate);
246 }
247
248 return maxComponents;
249 }
250
251 /** Performs Gradient optimization on a single atom using BarzilaiBorwein step width.
252 *
253 * \param _atom atom to anneal
254 * \param OldTimeStep old time step
255 * \param CurrentTimeStep current time step whose gradient we've just calculated
256 * \param TimeStepToSet time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
257 */
258 void annealAtom_BarzilaiBorwein(
259 atom * const _atom,
260 const int &OldTimeStep,
261 const int &CurrentTimeStep,
262 const int &TimeStepToSet
263 )
264 {
265 // atom's force vector gives steepest descent direction
266 const Vector &oldPosition = _atom->getPositionAtStep(OldTimeStep);
267 const Vector &currentPosition = _atom->getPositionAtStep(CurrentTimeStep);
268 const Vector &oldGradient = _atom->getAtomicForceAtStep(OldTimeStep);
269 const Vector &currentGradient = _atom->getAtomicForceAtStep(CurrentTimeStep);
270 LOG(4, "DEBUG: oldPosition for atom #" << _atom->getId() << " is " << oldPosition);
271 LOG(4, "DEBUG: currentPosition for atom #" << _atom->getId() << " is " << currentPosition);
272 LOG(4, "DEBUG: oldGradient for atom #" << _atom->getId() << " is " << oldGradient);
273 LOG(4, "DEBUG: currentGradient for atom #" << _atom->getId() << " is " << currentGradient);
274// LOG(4, "DEBUG: Force for atom #" << _atom->getId() << " is " << currentGradient);
275
276 // we use Barzilai-Borwein update with position reversed to get descent
277 const Vector PositionDifference = currentPosition - oldPosition;
278 const Vector GradientDifference = (currentGradient - oldGradient);
279 double stepwidth = getBarzilaiBorweinStepwidth(PositionDifference, GradientDifference);
280 Vector PositionUpdate = stepwidth * currentGradient;
281 LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
282
283 // finally set new values
284 _atom->setPositionAtStep(TimeStepToSet, currentPosition + PositionUpdate);
285 }
286
287 /** Performs Gradient optimization on the atoms using BarzilaiBorwein step width.
288 *
289 * \note this can only be called when there are at least two optimization
290 * time steps present, i.e. this must be preceded by a simple anneal().
291 *
292 * We assume that forces have just been calculated.
293 *
294 * \param _TimeStep time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
295 * \return to be filled with maximum force component over all atoms
296 */
297 Vector anneal_BarzilaiBorwein(
298 const int _TimeStep)
299 {
300 const int OldTimeStep = _TimeStep-2;
301 const int CurrentTimeStep = _TimeStep-1;
302 ASSERT( OldTimeStep >= 0,
303 "ForceAnnealing::anneal_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
304 ASSERT(currentStep > 1,
305 "ForceAnnealing::anneal_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
306
307 LOG(1, "STATUS: performing BarzilaiBorwein anneal at step #" << currentStep);
308
309 Vector maxComponents;
310 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
311 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
312
313 annealAtom_BarzilaiBorwein(*iter, OldTimeStep, CurrentTimeStep, _TimeStep);
314
315 // extract largest components for showing progress of annealing
316 const Vector &currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
317 for(size_t i=0;i<NDIM;++i)
318 maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
319 }
320
321 return maxComponents;
322 }
323
324 /** Performs Gradient optimization on the bonds with BarzilaiBorwein stepwdith.
325 *
326 * \note this can only be called when there are at least two optimization
327 * time steps present, i.e. this must be preceeded by a simple anneal().
328 *
329 * We assume that forces have just been calculated. These forces are projected
330 * onto the bonds and these are annealed subsequently by moving atoms in the
331 * bond neighborhood on either side conjunctively.
332 *
333 *
334 * \param _TimeStep time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
335 * \param maxComponents to be filled with maximum force component over all atoms
336 */
337 Vector annealWithBondGraph_BarzilaiBorwein(
338 const int _TimeStep)
339 {
340 const int OldTimeStep = _TimeStep-2;
341 const int CurrentTimeStep = _TimeStep-1;
342 ASSERT(OldTimeStep >= 0,
343 "annealWithBondGraph_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth, and the new one to update on already present.");
344 ASSERT(currentStep > 1,
345 "annealWithBondGraph_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
346
347 LOG(1, "STATUS: performing BarzilaiBorwein anneal on bonds at step #" << currentStep);
348
349 Vector maxComponents;
350
351 // get nodes on either side of selected bond via BFS discovery
352 BoostGraphCreator BGcreator;
353 BGcreator.createFromRange(
354 AtomicForceManipulator<T>::atoms.begin(),
355 AtomicForceManipulator<T>::atoms.end(),
356 AtomicForceManipulator<T>::atoms.size(),
357 BreadthFirstSearchGatherer::AlwaysTruePredicate);
358 BreadthFirstSearchGatherer NodeGatherer(BGcreator);
359
360 /** We assume that a force is local, i.e. a bond is too short yet and hence
361 * the atom needs to be moved. However, all the adjacent (bound) atoms might
362 * already be at the perfect distance. If we just move the atom alone, we ruin
363 * all the other bonds. Hence, it would be sensible to move every atom found
364 * through the bond graph in the direction of the force as well by the same
365 * PositionUpdate. This is almost what we are going to do, see below.
366 *
367 * This is to make the force a little more global in the sense of a multigrid
368 * solver that uses various coarser grids to transport errors more effectively
369 * over finely resolved grids.
370 *
371 */
372
373 /** The idea is that we project the gradients onto the bond vectors and determine
374 * from the sum of projected gradients from either side whether the bond is
375 * to contract or to expand. As the gradient acting as the normal vector of
376 * a plane supported at the position of the atom separates all bonds into two
377 * sets, we check whether all on one side are contracting and all on the other
378 * side are expanding. In this case we may move not only the atom itself but
379 * may propagate its update along a limited-horizon BFS to neighboring atoms.
380 *
381 */
382
383 // initialize helper class for bond vectors using bonds from range of atoms
384 BondVectors bv;
385 bv.setFromAtomRange< T >(
386 AtomicForceManipulator<T>::atoms.begin(),
387 AtomicForceManipulator<T>::atoms.end(),
388 _TimeStep); // use time step to update here as this is the current set of bonds
389
390 std::vector< // which bond side
391 std::vector<double> > // over all bonds
392 projected_forces; // one for leftatoms, one for rightatoms
393 projected_forces.resize(BondVectors::MAX_sides);
394 for (size_t j=0;j<BondVectors::MAX_sides;++j)
395 projected_forces[j].resize(bv.size(), 0.);
396
397 // for each atom we need to project the gradient
398 for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin();
399 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
400 const atom &walker = *(*iter);
401 const Vector &walkerGradient = walker.getAtomicForceAtStep(CurrentTimeStep);
402 const double GradientNorm = walkerGradient.Norm();
403 LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
404 << walker << " is " << walkerGradient << " with magnitude of "
405 << GradientNorm);
406
407 if (GradientNorm > MYEPSILON) {
408 bv.getProjectedGradientsForAtomAtStep(
409 walker, walkerGradient, CurrentTimeStep, projected_forces
410 );
411 } else {
412 LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
413 << MYEPSILON << " for atom " << walker);
414 // note that projected_forces is initialized to full length and filled
415 // with zeros. Hence, nothing to do here
416 }
417 }
418
419 std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
420 std::map<atomId_t, double> LargestUpdate_per_Atom; //!< check whether updates cancelled each other
421 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
422 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
423 atom &walker = *(*iter);
424
425 /// calculate step width
426 const Vector &oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
427 const Vector &currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
428 const Vector &oldGradient = (*iter)->getAtomicForceAtStep(OldTimeStep);
429 const Vector &currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
430 LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
431 LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition);
432 LOG(4, "DEBUG: oldGradient for atom #" << (*iter)->getId() << " is " << oldGradient);
433 LOG(4, "DEBUG: currentGradient for atom #" << (*iter)->getId() << " is " << currentGradient);
434// LOG(4, "DEBUG: Force for atom #" << (*iter)->getId() << " is " << currentGradient);
435
436 // we use Barzilai-Borwein update with position reversed to get descent
437 const Vector PositionDifference = currentPosition - oldPosition;
438 const Vector GradientDifference = (currentGradient - oldGradient);
439 double stepwidth = getBarzilaiBorweinStepwidth(PositionDifference, GradientDifference);
440 Vector PositionUpdate = stepwidth * currentGradient;
441 // cap updates (if non-zero) at 0.2 angstroem. BB tends to overshoot.
442 for (size_t i=0;i<NDIM;++i)
443 if (fabs(PositionUpdate[i]) > MYEPSILON)
444 PositionUpdate[i] = std::min(0.2, fabs(PositionUpdate[i]))*PositionUpdate[i]/fabs(PositionUpdate[i]);
445 LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
446
447 /** for each atom, we imagine a plane at the position of the atom with
448 * its atomic gradient as the normal vector. We go through all its bonds
449 * and check on which side of the plane the bond is. This defines whether
450 * the bond is contracting (+) or expanding (-) with respect to this atom.
451 *
452 * A bond has two atoms, however. Hence, we do this for either atom and
453 * look at the combination: Is it in sum contracting or expanding given
454 * both projected_forces?
455 */
456
457 /** go through all bonds and check projected_forces and side of plane
458 * the idea is that if all bonds on one side are contracting ones or expanding,
459 * respectively, then we may shift not only the atom with respect to its
460 * gradient but also its neighbors (towards contraction or towards
461 * expansion depending on direction of gradient).
462 * if they are mixed on both sides of the plane, then we simply shift
463 * only the atom itself.
464 * if they are not mixed on either side, then we also only shift the
465 * atom, namely away from expanding and towards contracting bonds.
466 *
467 * We may get this information right away by looking at the projected_forces.
468 * They give the atomic gradient of either atom projected onto the BondVector
469 * with an additional weight in [0,1].
470 */
471
472 // sign encodes side of plane and also encodes contracting(-) or expanding(+)
473 typedef std::vector<int> sides_t;
474 typedef std::vector<int> types_t;
475 sides_t sides;
476 types_t types;
477 const BondList& ListOfBonds = walker.getListOfBonds();
478 for(BondList::const_iterator bonditer = ListOfBonds.begin();
479 bonditer != ListOfBonds.end(); ++bonditer) {
480 const bond::ptr &current_bond = *bonditer;
481
482 // BondVector goes from bond::rightatom to bond::leftatom
483 const size_t index = bv.getIndexForBond(current_bond);
484 std::vector<double> &forcelist = (&walker == current_bond->leftatom) ?
485 projected_forces[BondVectors::leftside] : projected_forces[BondVectors::rightside];
486 // note that projected_forces has sign such as to indicate whether
487 // atomic gradient wants bond to contract (-) or expand (+).
488 // This goes into sides: Minus side points away from gradient, plus side point
489 // towards gradient.
490 //
491 // the sum of both bond sides goes into types, depending on which is
492 // stronger if either wants a different thing
493 const double &temp = forcelist[index];
494 if (fabs(temp) < MYEPSILON)
495 sides.push_back(1);
496 else
497 sides.push_back( -1.*temp/fabs(temp) ); // BondVectors has exactly opposite sign for sides decision
498 ASSERT( (sides.back() == 1) || (sides.back() == -1),
499 "ForceAnnealing() - sides is not in {-1,1}.");
500 const double sum =
501 projected_forces[BondVectors::leftside][index]+projected_forces[BondVectors::rightside][index];
502 types.push_back( sum/fabs(sum) );
503 LOG(4, "DEBUG: Bond " << *current_bond << " is on side " << sides.back()
504 << " and has type " << types.back());
505 }
506// /// check whether both conditions are compatible:
507// // i.e. either we have ++/-- for all entries in sides and types
508// // or we have +-/-+ for all entries
509// // hence, multiplying and taking the sum and its absolute value
510// // should be equal to the maximum number of entries
511// sides_t results;
512// std::transform(
513// sides.begin(), sides.end(),
514// types.begin(),
515// std::back_inserter(results),
516// std::multiplies<int>);
517// int result = abs(std::accumulate(results.begin(), results.end(), 0, std::plus<int>));
518
519 std::vector<size_t> first_per_side(2, (size_t)-1); //!< mark down one representative from either side
520 std::vector< std::vector<int> > types_per_side(2); //!< gather all types on each side
521 types_t::const_iterator typesiter = types.begin();
522 for (sides_t::const_iterator sidesiter = sides.begin();
523 sidesiter != sides.end(); ++sidesiter, ++typesiter) {
524 const size_t index = (*sidesiter+1)/2;
525 types_per_side[index].push_back(*typesiter);
526 if (first_per_side[index] == (size_t)-1)
527 first_per_side[index] = std::distance(const_cast<const sides_t &>(sides).begin(), sidesiter);
528 }
529 LOG(4, "DEBUG: First on side minus is " << first_per_side[0] << ", and first on side plus is "
530 << first_per_side[1]);
531 //!> enumerate types per side with a little witching with the numbers to allow easy setting from types
532 enum whichtypes_t {
533 contracting=0,
534 unset=1,
535 expanding=2,
536 mixed
537 };
538 std::vector<int> typeside(2, unset);
539 for(size_t i=0;i<2;++i) {
540 for (std::vector<int>::const_iterator tpsiter = types_per_side[i].begin();
541 tpsiter != types_per_side[i].end(); ++tpsiter) {
542 if (typeside[i] == unset) {
543 typeside[i] = *tpsiter+1; //contracting(0) or expanding(2)
544 } else {
545 if (typeside[i] != (*tpsiter+1)) // no longer he same type
546 typeside[i] = mixed;
547 }
548 }
549 }
550 LOG(4, "DEBUG: Minus side is " << typeside[0] << " and plus side is " << typeside[1]);
551
552 typedef std::vector< std::pair<atomId_t, atomId_t> > RemovedEdges_t;
553 if ((typeside[0] != mixed) || (typeside[1] != mixed)) {
554 const size_t sideno = ((typeside[0] != mixed) && (typeside[0] != unset)) ? 0 : 1;
555 LOG(4, "DEBUG: Chosen side is " << sideno << " with type " << typeside[sideno]);
556 ASSERT( (typeside[sideno] == contracting) || (typeside[sideno] == expanding),
557 "annealWithBondGraph_BB() - chosen side is neither expanding nor contracting.");
558 // one side is not mixed, all bonds on one side are of same type
559 // hence, find out which bonds to exclude
560 const BondList& ListOfBonds = walker.getListOfBonds();
561
562 // sideno is away (0) or in direction (1) of gradient
563 // tpyes[first_per_side[sideno]] is either contracting (-1) or expanding (+1)
564 // : side (i), where (i) means which bonds we keep for the BFS, bonds
565 // on side (-i) are removed
566 // If all bonds on side away (0) want expansion (+1), move towards side with atom: side 1
567 // if all bonds side towards (1) want contraction (-1), move away side with atom : side -1
568
569 // unsure whether this or do nothing in the remaining cases:
570 // If all bonds on side toward (1) want expansion (+1), move away side with atom : side -1
571 // (the reasoning is that the bond's other atom must have a stronger
572 // gradient in the same direction and they push along atoms in
573 // gradient direction: we don't want to interface with those.
574 // Hence, move atoms along on away side
575 // if all bonds side away (0) want contraction (-1), move towards side with atom: side 1
576 // (the reasoning is the same, don't interfere with update from
577 // stronger gradient)
578 // hence, the decision is only based on sides once we have picked a side
579 // depending on all bonds associated with have same good type.
580
581 // away from gradient (minus) and contracting
582 // or towards gradient (plus) and expanding
583 // gather all on same side and remove
584 const double sign =
585 (sides[first_per_side[sideno]] == types[first_per_side[sideno]])
586 ? sides[first_per_side[sideno]] : -1.*sides[first_per_side[sideno]];
587
588 LOG(4, "DEBUG: Removing edges from side with sign " << sign);
589 BondList::const_iterator bonditer = ListOfBonds.begin();
590 RemovedEdges_t RemovedEdges;
591 for (sides_t::const_iterator sidesiter = sides.begin();
592 sidesiter != sides.end(); ++sidesiter, ++bonditer) {
593 if (*sidesiter == sign) {
594 // remove the edge
595 const bond::ptr &current_bond = *bonditer;
596 LOG(5, "DEBUG: Removing edge " << *current_bond);
597 RemovedEdges.push_back( std::make_pair(
598 current_bond->leftatom->getId(),
599 current_bond->rightatom->getId())
600 );
601#ifndef NDEBUG
602 const bool status =
603#endif
604 BGcreator.removeEdge(RemovedEdges.back());
605 ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
606 }
607 }
608 // perform limited-horizon BFS
609 BoostGraphHelpers::Nodeset_t bondside_set;
610 BreadthFirstSearchGatherer::distance_map_t distance_map;
611 bondside_set = NodeGatherer(walker.getId(), max_distance);
612 distance_map = NodeGatherer.getDistances();
613 std::sort(bondside_set.begin(), bondside_set.end());
614
615 // re-add edge
616 for (RemovedEdges_t::const_iterator edgeiter = RemovedEdges.begin();
617 edgeiter != RemovedEdges.end(); ++edgeiter)
618 BGcreator.addEdge(edgeiter->first, edgeiter->second);
619
620 // update position with dampening factor on the discovered bonds
621 for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set.begin();
622 setiter != bondside_set.end(); ++setiter) {
623 const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
624 = distance_map.find(*setiter);
625 ASSERT( diter != distance_map.end(),
626 "ForceAnnealing() - could not find distance to an atom.");
627 const double factor = pow(damping_factor, diter->second+1);
628 LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
629 << factor << "*" << PositionUpdate);
630 if (GatheredUpdates.count((*setiter))) {
631 GatheredUpdates[(*setiter)] += factor*PositionUpdate;
632 LargestUpdate_per_Atom[(*setiter)] =
633 std::max(LargestUpdate_per_Atom[(*setiter)], factor*PositionUpdate.Norm());
634 } else {
635 GatheredUpdates.insert(
636 std::make_pair(
637 (*setiter),
638 factor*PositionUpdate) );
639 LargestUpdate_per_Atom.insert(
640 std::make_pair(
641 (*setiter),
642 factor*PositionUpdate.Norm()) );
643 }
644 }
645 } else {
646 // simple atomic annealing, i.e. damping factor of 1
647 LOG(3, "DEBUG: Update for atom #" << walker.getId() << " will be " << PositionUpdate);
648 GatheredUpdates.insert(
649 std::make_pair(
650 walker.getId(),
651 PositionUpdate) );
652 LargestUpdate_per_Atom.insert(
653 std::make_pair(
654 walker.getId(),
655 PositionUpdate.Norm()) );
656 }
657 }
658
659 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
660 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
661 atom &walker = *(*iter);
662 // extract largest components for showing progress of annealing
663 const Vector &currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep);
664 for(size_t i=0;i<NDIM;++i)
665 maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
666 }
667
668// // remove center of weight translation from gathered updates
669// Vector CommonTranslation;
670// for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
671// iter != GatheredUpdates.end(); ++iter) {
672// const Vector &update = iter->second;
673// CommonTranslation += update;
674// }
675// CommonTranslation *= 1./(double)GatheredUpdates.size();
676// LOG(3, "DEBUG: Subtracting common translation " << CommonTranslation
677// << " from all updates.");
678
679 // apply the gathered updates and set remnant gradients for atomic annealing
680 Vector LargestUpdate;
681 for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
682 iter != GatheredUpdates.end(); ++iter) {
683 const atomId_t &atomid = iter->first;
684 const Vector &update = iter->second;
685 atom* const walker = World::getInstance().getAtom(AtomById(atomid));
686 ASSERT( walker != NULL,
687 "ForceAnnealing() - walker with id "+toString(atomid)+" has suddenly disappeared.");
688 LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
689 << ", namely " << *walker);
690 for (size_t i=0;i<NDIM;++i)
691 LargestUpdate[i] = std::max(LargestUpdate[i], fabs(update[i]));
692
693 std::map<atomId_t, double>::const_iterator largestiter = LargestUpdate_per_Atom.find(atomid);
694 ASSERT( largestiter != LargestUpdate_per_Atom.end(),
695 "ForceAnnealing() - walker with id "+toString(atomid)+" not in LargestUpdates.");
696 // if we had large updates but their sum is very small
697 if (update.Norm()/largestiter->second > MYEPSILON) {
698 walker->setPositionAtStep(_TimeStep,
699 walker->getPositionAtStep(CurrentTimeStep) + update); // - CommonTranslation);
700 } else {
701 // then recalc update with simple anneal
702 LOG(2, "WARNING: Updates on atom " << *iter << " cancel themselves, performing simple anneal step.");
703 annealAtom_BarzilaiBorwein(walker, OldTimeStep, CurrentTimeStep, _TimeStep);
704 }
705 }
706 LOG(1, "STATUS: Largest absolute update components are " << LargestUpdate);
707
708 return maxComponents;
709 }
710
711 /** Reset function to unset static entities and artificial velocities.
712 *
713 */
714 void reset()
715 {
716 currentDeltat = 0.;
717 currentStep = 0;
718 }
719
720private:
721 //!> contains the current step in relation to maxsteps
722 static size_t currentStep;
723 //!> contains the maximum number of steps, determines initial and final step with currentStep
724 size_t maxSteps;
725 static double currentDeltat;
726 //!> minimum deltat for internal while loop (adaptive step width)
727 static double MinimumDeltat;
728 //!> contains the maximum bond graph distance up to which shifts of a single atom are spread
729 const int max_distance;
730 //!> the shifted is dampened by this factor with the power of the bond graph distance to the shift causing atom
731 const double damping_factor;
732 //!> threshold for force components to stop annealing
733 const double FORCE_THRESHOLD;
734};
735
736template <class T>
737double ForceAnnealing<T>::currentDeltat = 0.;
738template <class T>
739size_t ForceAnnealing<T>::currentStep = 0;
740template <class T>
741double ForceAnnealing<T>::MinimumDeltat = 1e-8;
742
743#endif /* FORCEANNEALING_HPP_ */
Note: See TracBrowser for help on using the repository browser.