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 | */
|
---|
38 | template <class T>
|
---|
39 | class ForceAnnealing : public AtomicForceManipulator<T>
|
---|
40 | {
|
---|
41 | public:
|
---|
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 |
|
---|
165 | private:
|
---|
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 |
|
---|
175 | template <class T>
|
---|
176 | double ForceAnnealing<T>::currentDeltat = 0.;
|
---|
177 | template <class T>
|
---|
178 | size_t ForceAnnealing<T>::currentStep = 0;
|
---|
179 | template <class T>
|
---|
180 | double ForceAnnealing<T>::MinimumDeltat = 1e-8;
|
---|
181 |
|
---|
182 | #endif /* FORCEANNEALING_HPP_ */
|
---|