source: src/Dynamics/ForceAnnealing.hpp@ 322d58

AutomationFragmentation_failures Candidate_v1.6.1 ChemicalSpaceEvaluator Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Exclude_Hydrogens_annealWithBondGraph ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_contraction-expansion Gui_displays_atomic_force_velocity JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool PythonUI_with_named_parameters StoppableMakroAction TremoloParser_IncreasedPrecision
Last change on this file since 322d58 was 322d58, checked in by Frederik Heber <frederik.heber@…>, 7 years ago

ForceAnnealing now uses step width according to Barzilai-Borwein method.

  • this is not as good as what MPQC does but as long as Actions don't have an internal state, we cannot do anything better, e.g. a real line-search.
  • removed paramater delta from optimize-structure.
  • DOCU: Updated entry in user guide and noted Barzilai-Borwein method.
  • TESTS: removed deltat from force annealing regression tests.
  • TESTS: switched off StructureOptimization integration tests as poolworkers crash at the moment, making the tests hang indefinitely.
  • Property mode set to 100644
File size: 6.5 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 "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 "Dynamics/AtomicForceManipulator.hpp"
23#include "Fragmentation/ForceMatrix.hpp"
24#include "Helpers/helpers.hpp"
25#include "Helpers/defs.hpp"
26#include "LinearAlgebra/Vector.hpp"
27#include "Thermostats/ThermoStatContainer.hpp"
28#include "Thermostats/Thermostat.hpp"
29#include "World.hpp"
30
31/** This class is the essential build block for performing structual optimization.
32 *
33 * Sadly, we have to use some static instances as so far values cannot be passed
34 * between actions. Hence, we need to store the current step and the adaptive-
35 * step width (we cannot perform a linesearch, as we have no control over the
36 * calculation of the forces).
37 */
38template <class T>
39class ForceAnnealing : public AtomicForceManipulator<T>
40{
41public:
42 /** Constructor of class ForceAnnealing.
43 *
44 * \note We use a fixed delta t of 1.
45 *
46 * \param _atoms set of atoms to integrate
47 * \param _Deltat time step width in atomic units
48 * \param _IsAngstroem whether length units are in angstroem or bohr radii
49 * \param _maxSteps number of optimization steps to perform
50 */
51 ForceAnnealing(
52 AtomSetMixin<T> &_atoms,
53 bool _IsAngstroem,
54 const size_t _maxSteps) :
55 AtomicForceManipulator<T>(_atoms, 1., _IsAngstroem),
56 maxSteps(_maxSteps)
57 {}
58 /** Destructor of class ForceAnnealing.
59 *
60 */
61 ~ForceAnnealing()
62 {}
63
64 /** Performs Gradient optimization.
65 *
66 * We assume that forces have just been calculated.
67 *
68 *
69 * \param NextStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
70 * \param offset offset in matrix file to the first force component
71 * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
72 */
73 void operator()(const int NextStep, const size_t offset)
74 {
75 // make sum of forces equal zero
76 AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(offset,NextStep);
77
78 // are we in initial step? Then set static entities
79 if (currentStep == 0) {
80 currentDeltat = AtomicForceManipulator<T>::Deltat;
81 currentStep = 1;
82 LOG(2, "DEBUG: Initial step, setting values, current step is #" << currentStep);
83 } else {
84 ++currentStep;
85 LOG(2, "DEBUG: current step is #" << currentStep);
86 }
87
88 Vector maxComponents(zeroVec);
89 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
90 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
91 // atom's force vector gives steepest descent direction
92 const Vector oldPosition = (*iter)->getPositionAtStep(NextStep-2 >= 0 ? NextStep - 2 : 0);
93 const Vector currentPosition = (*iter)->getPosition();
94 const Vector oldGradient = (*iter)->getAtomicForceAtStep(NextStep-2 >= 0 ? NextStep - 2 : 0);
95 const Vector currentGradient = (*iter)->getAtomicForce();
96 LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
97
98 // we use Barzilai-Borwein update with position reversed to get descent
99 const Vector GradientDifference = (currentGradient - oldGradient);
100 const double stepwidth =
101 fabs((currentPosition - oldPosition).ScalarProduct(GradientDifference))/
102 GradientDifference.NormSquared();
103 Vector PositionUpdate = stepwidth * currentGradient;
104 if (fabs(stepwidth) < 1e-10) {
105 // dont' warn in first step, deltat usage normal
106 if (currentStep != 1)
107 ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
108 PositionUpdate = currentDeltat * currentGradient;
109 }
110 LOG(3, "DEBUG: Update would be " << PositionUpdate);
111
112 // extract largest components for showing progress of annealing
113 for(size_t i=0;i<NDIM;++i)
114 if (currentGradient[i] > maxComponents[i])
115 maxComponents[i] = currentGradient[i];
116
117 // are we in initial step? Then don't check against velocity
118 if ((currentStep > 1) && (!(*iter)->getAtomicVelocity().IsZero()))
119 // update with currentDelta tells us how the current gradient relates to
120 // the last one: If it has become larger, reduce currentDelta
121 if ((PositionUpdate.ScalarProduct((*iter)->getAtomicVelocity()) < 0)
122 && (currentDeltat > MinimumDeltat)) {
123 currentDeltat = .5*currentDeltat;
124 LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate.NormSquared()
125 << " > " << (*iter)->getAtomicVelocity().NormSquared()
126 << ", decreasing deltat: " << currentDeltat);
127 PositionUpdate = currentDeltat * currentGradient;
128 }
129
130 // finally set new values
131 (*iter)->setPosition(currentPosition + PositionUpdate);
132 (*iter)->setAtomicVelocity(PositionUpdate);
133 //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
134// (*iter)->VelocityVerletUpdateU((*iter)->getId(), NextStep-1, Deltat, IsAngstroem);
135
136 // reset force vector for next step except on final one
137 if (currentStep != maxSteps)
138 (*iter)->setAtomicForce(zeroVec);
139 }
140
141 LOG(1, "STATUS: Largest remaining force components at step #"
142 << currentStep << " are " << maxComponents);
143
144 // are we in final step? Remember to reset static entities
145 if (currentStep == maxSteps) {
146 LOG(2, "DEBUG: Final step, resetting values");
147 reset();
148 }
149 }
150
151 /** Reset function to unset static entities and artificial velocities.
152 *
153 */
154 void reset()
155 {
156 currentDeltat = 0.;
157 currentStep = 0;
158
159 // reset (artifical) velocities
160 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
161 iter != AtomicForceManipulator<T>::atoms.end(); ++iter)
162 (*iter)->setAtomicVelocity(zeroVec);
163 }
164
165private:
166 //!> contains the current step in relation to maxsteps
167 static size_t currentStep;
168 //!> contains the maximum number of steps, determines initial and final step with currentStep
169 size_t maxSteps;
170 static double currentDeltat;
171 //!> minimum deltat for internal while loop (adaptive step width)
172 static double MinimumDeltat;
173};
174
175template <class T>
176double ForceAnnealing<T>::currentDeltat = 0.;
177template <class T>
178size_t ForceAnnealing<T>::currentStep = 0;
179template <class T>
180double ForceAnnealing<T>::MinimumDeltat = 1e-8;
181
182#endif /* FORCEANNEALING_HPP_ */
Note: See TracBrowser for help on using the repository browser.