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 |
|
---|
20 | #include <boost/bind.hpp>
|
---|
21 |
|
---|
22 | #include "Atom/atom.hpp"
|
---|
23 | #include "Atom/AtomSet.hpp"
|
---|
24 | #include "CodePatterns/Assert.hpp"
|
---|
25 | #include "CodePatterns/Info.hpp"
|
---|
26 | #include "CodePatterns/Log.hpp"
|
---|
27 | #include "CodePatterns/Verbose.hpp"
|
---|
28 | #include "Descriptors/AtomIdDescriptor.hpp"
|
---|
29 | #include "Dynamics/AtomicForceManipulator.hpp"
|
---|
30 | #include "Dynamics/BondVectors.hpp"
|
---|
31 | #include "Fragmentation/ForceMatrix.hpp"
|
---|
32 | #include "Graph/BoostGraphCreator.hpp"
|
---|
33 | #include "Graph/BoostGraphHelpers.hpp"
|
---|
34 | #include "Graph/BreadthFirstSearchGatherer.hpp"
|
---|
35 | #include "Helpers/helpers.hpp"
|
---|
36 | #include "Helpers/defs.hpp"
|
---|
37 | #include "LinearAlgebra/LinearSystemOfEquations.hpp"
|
---|
38 | #include "LinearAlgebra/MatrixContent.hpp"
|
---|
39 | #include "LinearAlgebra/Vector.hpp"
|
---|
40 | #include "LinearAlgebra/VectorContent.hpp"
|
---|
41 | #include "Thermostats/ThermoStatContainer.hpp"
|
---|
42 | #include "Thermostats/Thermostat.hpp"
|
---|
43 | #include "World.hpp"
|
---|
44 |
|
---|
45 | /** This class is the essential build block for performing structural optimization.
|
---|
46 | *
|
---|
47 | * Sadly, we have to use some static instances as so far values cannot be passed
|
---|
48 | * between actions. Hence, we need to store the current step and the adaptive-
|
---|
49 | * step width (we cannot perform a line search, as we have no control over the
|
---|
50 | * calculation of the forces).
|
---|
51 | *
|
---|
52 | * However, we do use the bond graph, i.e. if a single atom needs to be shifted
|
---|
53 | * to the left, then the whole molecule left of it is shifted, too. This is
|
---|
54 | * controlled by the \a max_distance parameter.
|
---|
55 | */
|
---|
56 | template <class T>
|
---|
57 | class ForceAnnealing : public AtomicForceManipulator<T>
|
---|
58 | {
|
---|
59 | public:
|
---|
60 | /** Constructor of class ForceAnnealing.
|
---|
61 | *
|
---|
62 | * \note We use a fixed delta t of 1.
|
---|
63 | *
|
---|
64 | * \param _atoms set of atoms to integrate
|
---|
65 | * \param _Deltat time step width in atomic units
|
---|
66 | * \param _IsAngstroem whether length units are in angstroem or bohr radii
|
---|
67 | * \param _maxSteps number of optimization steps to perform
|
---|
68 | * \param _max_distance up to this bond order is bond graph taken into account.
|
---|
69 | */
|
---|
70 | ForceAnnealing(
|
---|
71 | AtomSetMixin<T> &_atoms,
|
---|
72 | const double _Deltat,
|
---|
73 | bool _IsAngstroem,
|
---|
74 | const size_t _maxSteps,
|
---|
75 | const int _max_distance,
|
---|
76 | const double _damping_factor) :
|
---|
77 | AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem),
|
---|
78 | maxSteps(_maxSteps),
|
---|
79 | max_distance(_max_distance),
|
---|
80 | damping_factor(_damping_factor)
|
---|
81 | {}
|
---|
82 |
|
---|
83 | /** Destructor of class ForceAnnealing.
|
---|
84 | *
|
---|
85 | */
|
---|
86 | ~ForceAnnealing()
|
---|
87 | {}
|
---|
88 |
|
---|
89 | /** Performs Gradient optimization.
|
---|
90 | *
|
---|
91 | * We assume that forces have just been calculated.
|
---|
92 | *
|
---|
93 | *
|
---|
94 | * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
|
---|
95 | * \param offset offset in matrix file to the first force component
|
---|
96 | * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
|
---|
97 | */
|
---|
98 | void operator()(
|
---|
99 | const int _CurrentTimeStep,
|
---|
100 | const size_t _offset,
|
---|
101 | const bool _UseBondgraph)
|
---|
102 | {
|
---|
103 | // make sum of forces equal zero
|
---|
104 | AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(
|
---|
105 | _offset,
|
---|
106 | _CurrentTimeStep-1>=0 ? _CurrentTimeStep - 1 : 0);
|
---|
107 |
|
---|
108 | // are we in initial step? Then set static entities
|
---|
109 | Vector maxComponents(zeroVec);
|
---|
110 | if (currentStep == 0) {
|
---|
111 | currentDeltat = AtomicForceManipulator<T>::Deltat;
|
---|
112 | currentStep = 1;
|
---|
113 | LOG(2, "DEBUG: Initial step, setting values, current step is #" << currentStep);
|
---|
114 |
|
---|
115 | // always use atomic annealing on first step
|
---|
116 | anneal(_CurrentTimeStep, _offset, maxComponents);
|
---|
117 | } else {
|
---|
118 | ++currentStep;
|
---|
119 | LOG(2, "DEBUG: current step is #" << currentStep);
|
---|
120 |
|
---|
121 | if (_UseBondgraph)
|
---|
122 | annealWithBondGraph(_CurrentTimeStep, _offset, maxComponents);
|
---|
123 | else
|
---|
124 | anneal(_CurrentTimeStep, _offset, maxComponents);
|
---|
125 | }
|
---|
126 |
|
---|
127 |
|
---|
128 | LOG(1, "STATUS: Largest remaining force components at step #"
|
---|
129 | << currentStep << " are " << maxComponents);
|
---|
130 |
|
---|
131 | // are we in final step? Remember to reset static entities
|
---|
132 | if (currentStep == maxSteps) {
|
---|
133 | LOG(2, "DEBUG: Final step, resetting values");
|
---|
134 | reset();
|
---|
135 | }
|
---|
136 | }
|
---|
137 |
|
---|
138 | /** Helper function to calculate the Barzilai-Borwein stepwidth.
|
---|
139 | *
|
---|
140 | * \param _PositionDifference difference in position between current and last step
|
---|
141 | * \param _GradientDifference difference in gradient between current and last step
|
---|
142 | * \return step width according to Barzilai-Borwein
|
---|
143 | */
|
---|
144 | double getBarzilaiBorweinStepwidth(const Vector &_PositionDifference, const Vector &_GradientDifference)
|
---|
145 | {
|
---|
146 | double stepwidth = 0.;
|
---|
147 | if (_GradientDifference.NormSquared() > MYEPSILON)
|
---|
148 | stepwidth = fabs(_PositionDifference.ScalarProduct(_GradientDifference))/
|
---|
149 | _GradientDifference.NormSquared();
|
---|
150 | if (fabs(stepwidth) < 1e-10) {
|
---|
151 | // dont' warn in first step, deltat usage normal
|
---|
152 | if (currentStep != 1)
|
---|
153 | ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
|
---|
154 | stepwidth = currentDeltat;
|
---|
155 | }
|
---|
156 | return stepwidth;
|
---|
157 | }
|
---|
158 |
|
---|
159 | /** Performs Gradient optimization on the atoms.
|
---|
160 | *
|
---|
161 | * We assume that forces have just been calculated.
|
---|
162 | *
|
---|
163 | * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
|
---|
164 | * \param offset offset in matrix file to the first force component
|
---|
165 | * \param maxComponents to be filled with maximum force component over all atoms
|
---|
166 | */
|
---|
167 | void anneal(
|
---|
168 | const int CurrentTimeStep,
|
---|
169 | const size_t offset,
|
---|
170 | Vector &maxComponents)
|
---|
171 | {
|
---|
172 | bool deltat_decreased = false;
|
---|
173 | for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
|
---|
174 | iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
|
---|
175 | // atom's force vector gives steepest descent direction
|
---|
176 | const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
|
---|
177 | const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
|
---|
178 | const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
|
---|
179 | const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
|
---|
180 | LOG(4, "DEBUG: oldPosition for atom " << **iter << " is " << oldPosition);
|
---|
181 | LOG(4, "DEBUG: currentPosition for atom " << **iter << " is " << currentPosition);
|
---|
182 | LOG(4, "DEBUG: oldGradient for atom " << **iter << " is " << oldGradient);
|
---|
183 | LOG(4, "DEBUG: currentGradient for atom " << **iter << " is " << currentGradient);
|
---|
184 | // LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
|
---|
185 |
|
---|
186 | // we use Barzilai-Borwein update with position reversed to get descent
|
---|
187 | const double stepwidth = getBarzilaiBorweinStepwidth(
|
---|
188 | currentPosition - oldPosition, currentGradient - oldGradient);
|
---|
189 | Vector PositionUpdate = stepwidth * currentGradient;
|
---|
190 | LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
|
---|
191 |
|
---|
192 | // extract largest components for showing progress of annealing
|
---|
193 | for(size_t i=0;i<NDIM;++i)
|
---|
194 | maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
|
---|
195 |
|
---|
196 | // steps may go back and forth again (updates are of same magnitude but
|
---|
197 | // have different sign: Check whether this is the case and one step with
|
---|
198 | // deltat to interrupt this sequence
|
---|
199 | const Vector PositionDifference = currentPosition - oldPosition;
|
---|
200 | if ((currentStep > 1) && (!PositionDifference.IsZero()))
|
---|
201 | if ((PositionUpdate.ScalarProduct(PositionDifference) < 0)
|
---|
202 | && (fabs(PositionUpdate.NormSquared()-PositionDifference.NormSquared()) < 1e-3)) {
|
---|
203 | // for convergence we want a null sequence here, too
|
---|
204 | if (!deltat_decreased) {
|
---|
205 | deltat_decreased = true;
|
---|
206 | currentDeltat = .5*currentDeltat;
|
---|
207 | }
|
---|
208 | LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate
|
---|
209 | << " > " << PositionDifference
|
---|
210 | << ", using deltat: " << currentDeltat);
|
---|
211 | PositionUpdate = currentDeltat * currentGradient;
|
---|
212 | }
|
---|
213 |
|
---|
214 | // finally set new values
|
---|
215 | (*iter)->setPosition(currentPosition + PositionUpdate);
|
---|
216 | }
|
---|
217 | }
|
---|
218 |
|
---|
219 | /** Performs Gradient optimization on the bonds.
|
---|
220 | *
|
---|
221 | * We assume that forces have just been calculated. These forces are projected
|
---|
222 | * onto the bonds and these are annealed subsequently by moving atoms in the
|
---|
223 | * bond neighborhood on either side conjunctively.
|
---|
224 | *
|
---|
225 | *
|
---|
226 | * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
|
---|
227 | * \param offset offset in matrix file to the first force component
|
---|
228 | * \param maxComponents to be filled with maximum force component over all atoms
|
---|
229 | */
|
---|
230 | void annealWithBondGraph(
|
---|
231 | const int CurrentTimeStep,
|
---|
232 | const size_t offset,
|
---|
233 | Vector &maxComponents)
|
---|
234 | {
|
---|
235 | // get nodes on either side of selected bond via BFS discovery
|
---|
236 | BoostGraphCreator BGcreator;
|
---|
237 | BGcreator.createFromRange(
|
---|
238 | AtomicForceManipulator<T>::atoms.begin(),
|
---|
239 | AtomicForceManipulator<T>::atoms.end(),
|
---|
240 | AtomicForceManipulator<T>::atoms.size(),
|
---|
241 | BreadthFirstSearchGatherer::AlwaysTruePredicate);
|
---|
242 | BreadthFirstSearchGatherer NodeGatherer(BGcreator);
|
---|
243 |
|
---|
244 | /// We assume that a force is local, i.e. a bond is too short yet and hence
|
---|
245 | /// the atom needs to be moved. However, all the adjacent (bound) atoms might
|
---|
246 | /// already be at the perfect distance. If we just move the atom alone, we ruin
|
---|
247 | /// all the other bonds. Hence, it would be sensible to move every atom found
|
---|
248 | /// through the bond graph in the direction of the force as well by the same
|
---|
249 | /// PositionUpdate. This is almost what we are going to do.
|
---|
250 |
|
---|
251 | /// One issue is: If we need to shorten bond, then we use the PositionUpdate
|
---|
252 | /// also on the the other bond partner already. This is because it is in the
|
---|
253 | /// direction of the bond. Therefore, the update is actually performed twice on
|
---|
254 | /// each bond partner, i.e. the step size is twice as large as it should be.
|
---|
255 | /// This problem only occurs when bonds need to be shortened, not when they
|
---|
256 | /// need to be made longer (then the force vector is facing the other
|
---|
257 | /// direction than the bond vector).
|
---|
258 | /// As a remedy we need to average the force on either end of the bond and
|
---|
259 | /// check whether each gradient points inwards out or outwards with respect
|
---|
260 | /// to the bond and then shift accordingly.
|
---|
261 |
|
---|
262 | /// One more issue is that the projection onto the bond directions does not
|
---|
263 | /// recover the gradient but may be larger as the bond directions are a
|
---|
264 | /// generating system and not a basis (e.g. 3 bonds on a plane where 2 would
|
---|
265 | /// suffice to span the plane). To this end, we need to account for the
|
---|
266 | /// overestimation and obtain a weighting for each bond.
|
---|
267 |
|
---|
268 | // initialize helper class for bond vectors using bonds from range of atoms
|
---|
269 | BondVectors bv;
|
---|
270 | bv.setFromAtomRange< T >(
|
---|
271 | AtomicForceManipulator<T>::atoms.begin(),
|
---|
272 | AtomicForceManipulator<T>::atoms.end(),
|
---|
273 | CurrentTimeStep);
|
---|
274 | const BondVectors::container_t &sorted_bonds = bv.getSorted();
|
---|
275 |
|
---|
276 | // knowing the number of bonds in total, we can setup the storage for the
|
---|
277 | // projected forces
|
---|
278 | enum whichatom_t {
|
---|
279 | leftside=0,
|
---|
280 | rightside=1,
|
---|
281 | MAX_sides
|
---|
282 | };
|
---|
283 | std::vector< // time step
|
---|
284 | std::vector< // which bond side
|
---|
285 | std::vector<double> > // over all bonds
|
---|
286 | > projected_forces(2); // one for leftatoms, one for rightatoms (and for both time steps)
|
---|
287 | for (size_t i=0;i<2;++i) {
|
---|
288 | projected_forces[i].resize(MAX_sides);
|
---|
289 | for (size_t j=0;j<MAX_sides;++j)
|
---|
290 | projected_forces[i][j].resize(sorted_bonds.size(), 0.);
|
---|
291 | }
|
---|
292 |
|
---|
293 | // for each atom we need to gather weights and then project the gradient
|
---|
294 | typedef std::deque<double> weights_t;
|
---|
295 | typedef std::map<atomId_t, weights_t > weights_per_atom_t;
|
---|
296 | std::vector<weights_per_atom_t> weights_per_atom(2);
|
---|
297 | for (size_t timestep = 0; timestep <= 1; ++timestep) {
|
---|
298 | const size_t CurrentStep = CurrentTimeStep-timestep-1 >= 0 ? CurrentTimeStep-timestep-1 : 0;
|
---|
299 | LOG(2, "DEBUG: CurrentTimeStep is " << CurrentTimeStep
|
---|
300 | << ", timestep is " << timestep
|
---|
301 | << ", and CurrentStep is " << CurrentStep);
|
---|
302 |
|
---|
303 | // get all bond vectors for this time step (from the perspective of the
|
---|
304 | // bonds taken from the currentStep)
|
---|
305 | const BondVectors::mapped_t bondvectors = bv.getBondVectorsAtStep(CurrentStep);
|
---|
306 |
|
---|
307 | for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin();
|
---|
308 | iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
|
---|
309 | const atom &walker = *(*iter);
|
---|
310 | const Vector &walkerGradient = walker.getAtomicForceAtStep(CurrentStep);
|
---|
311 | LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
|
---|
312 | << walker << " is " << walkerGradient << " with magnitude of "
|
---|
313 | << walkerGradient.Norm());
|
---|
314 |
|
---|
315 | const BondList& ListOfBonds = walker.getListOfBonds();
|
---|
316 | if (walkerGradient.Norm() > MYEPSILON) {
|
---|
317 |
|
---|
318 | // gather subset of BondVectors for the current atom
|
---|
319 | std::vector<Vector> BondVectors;
|
---|
320 | for(BondList::const_iterator bonditer = ListOfBonds.begin();
|
---|
321 | bonditer != ListOfBonds.end(); ++bonditer) {
|
---|
322 | const bond::ptr ¤t_bond = *bonditer;
|
---|
323 | const BondVectors::mapped_t::const_iterator bviter =
|
---|
324 | bondvectors.find(current_bond);
|
---|
325 | ASSERT( bviter != bondvectors.end(),
|
---|
326 | "ForceAnnealing() - cannot find current_bond ?");
|
---|
327 | ASSERT( bviter != bondvectors.end(),
|
---|
328 | "ForceAnnealing - cannot find current bond "+toString(*current_bond)
|
---|
329 | +" in bonds.");
|
---|
330 | BondVectors.push_back(bviter->second);
|
---|
331 | }
|
---|
332 | LOG(4, "DEBUG: BondVectors for atom #" << walker.getId() << ": " << BondVectors);
|
---|
333 |
|
---|
334 | // go through all its bonds and calculate what magnitude is represented
|
---|
335 | // by the others i.e. sum of scalar products against other bonds
|
---|
336 | std::pair<weights_per_atom_t::iterator, bool> inserter =
|
---|
337 | weights_per_atom[timestep].insert(
|
---|
338 | std::make_pair(walker.getId(), weights_t()) );
|
---|
339 | ASSERT( inserter.second,
|
---|
340 | "ForceAnnealing::operator() - weight map for atom "+toString(walker)
|
---|
341 | +" and time step "+toString(timestep)+" already filled?");
|
---|
342 | weights_t &weights = inserter.first->second;
|
---|
343 | for (std::vector<Vector>::const_iterator iter = BondVectors.begin();
|
---|
344 | iter != BondVectors.end(); ++iter) {
|
---|
345 | std::vector<double> scps;
|
---|
346 | scps.reserve(BondVectors.size());
|
---|
347 | std::transform(
|
---|
348 | BondVectors.begin(), BondVectors.end(),
|
---|
349 | std::back_inserter(scps),
|
---|
350 | boost::bind(static_cast< double (*)(double) >(&fabs),
|
---|
351 | boost::bind(&Vector::ScalarProduct, boost::cref(*iter), _1))
|
---|
352 | );
|
---|
353 | const double scp_sum = std::accumulate(scps.begin(), scps.end(), 0.);
|
---|
354 | ASSERT( (scp_sum-1.) > -MYEPSILON,
|
---|
355 | "ForceAnnealing() - sum of weights must be equal or larger one but is "
|
---|
356 | +toString(scp_sum));
|
---|
357 | weights.push_back( 1./scp_sum );
|
---|
358 | }
|
---|
359 | LOG(4, "DEBUG: Weights for atom #" << walker.getId() << ": " << weights);
|
---|
360 |
|
---|
361 | // for testing we check whether all weighted scalar products now come out as 1.
|
---|
362 | #ifndef NDEBUG
|
---|
363 | for (std::vector<Vector>::const_iterator iter = BondVectors.begin();
|
---|
364 | iter != BondVectors.end(); ++iter) {
|
---|
365 | std::vector<double> scps;
|
---|
366 | scps.reserve(BondVectors.size());
|
---|
367 | std::transform(
|
---|
368 | BondVectors.begin(), BondVectors.end(),
|
---|
369 | weights.begin(),
|
---|
370 | std::back_inserter(scps),
|
---|
371 | boost::bind(static_cast< double (*)(double) >(&fabs),
|
---|
372 | boost::bind(std::multiplies<double>(),
|
---|
373 | boost::bind(&Vector::ScalarProduct, boost::cref(*iter), _1),
|
---|
374 | _2))
|
---|
375 | );
|
---|
376 | const double scp_sum = std::accumulate(scps.begin(), scps.end(), 0.);
|
---|
377 | ASSERT( fabs(scp_sum - 1.) < MYEPSILON,
|
---|
378 | "ForceAnnealing::operator() - for BondVector "+toString(*iter)
|
---|
379 | +" we have weighted scalar product of "+toString(scp_sum)+" != 1.");
|
---|
380 | }
|
---|
381 | #endif
|
---|
382 |
|
---|
383 | // projected gradient over all bonds and place in one of projected_forces
|
---|
384 | // using the obtained weights
|
---|
385 | {
|
---|
386 | weights_t::const_iterator weightiter = weights.begin();
|
---|
387 | std::vector<Vector>::const_iterator vectoriter = BondVectors.begin();
|
---|
388 | Vector forcesum_check;
|
---|
389 | for(BondList::const_iterator bonditer = ListOfBonds.begin();
|
---|
390 | bonditer != ListOfBonds.end(); ++bonditer, ++weightiter, ++vectoriter) {
|
---|
391 | const bond::ptr ¤t_bond = *bonditer;
|
---|
392 | const Vector &BondVector = *vectoriter;
|
---|
393 |
|
---|
394 | std::vector<double> &forcelist = (&walker == current_bond->leftatom) ?
|
---|
395 | projected_forces[timestep][leftside] : projected_forces[timestep][rightside];
|
---|
396 | const size_t index = bv.getIndexForBond(current_bond);
|
---|
397 | ASSERT( index != (size_t)-1,
|
---|
398 | "ForceAnnealing() - could not find bond "+toString(*current_bond)
|
---|
399 | +" in bondvectors");
|
---|
400 | forcelist[index] = (*weightiter)*walkerGradient.ScalarProduct(BondVector);
|
---|
401 | LOG(4, "DEBUG: BondVector " << BondVector << " receives projected force of "
|
---|
402 | << forcelist[index]);
|
---|
403 | forcesum_check += forcelist[index] * BondVector;
|
---|
404 | }
|
---|
405 | ASSERT( weightiter == weights.end(),
|
---|
406 | "ForceAnnealing - weightiter is not at end when it should be.");
|
---|
407 | ASSERT( vectoriter == BondVectors.end(),
|
---|
408 | "ForceAnnealing - vectoriter is not at end when it should be.");
|
---|
409 | LOG(3, "DEBUG: sum of projected forces is " << forcesum_check);
|
---|
410 | }
|
---|
411 |
|
---|
412 | } else {
|
---|
413 | LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
|
---|
414 | << MYEPSILON << " for atom " << walker);
|
---|
415 | // note that projected_forces is initialized to full length and filled
|
---|
416 | // with zeros. Hence, nothing to do here
|
---|
417 | }
|
---|
418 | }
|
---|
419 | }
|
---|
420 |
|
---|
421 | // step through each bond and shift the atoms
|
---|
422 | std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
|
---|
423 |
|
---|
424 | LOG(3, "DEBUG: current step is " << currentStep << ", given time step is " << CurrentTimeStep);
|
---|
425 | const BondVectors::mapped_t bondvectors = bv.getBondVectorsAtStep(CurrentTimeStep);
|
---|
426 |
|
---|
427 | for (BondVectors::container_t::const_iterator bondsiter = sorted_bonds.begin();
|
---|
428 | bondsiter != sorted_bonds.end(); ++bondsiter) {
|
---|
429 | const bond::ptr ¤t_bond = *bondsiter;
|
---|
430 | const size_t index = bv.getIndexForBond(current_bond);
|
---|
431 | const atom* bondatom[MAX_sides] = {
|
---|
432 | current_bond->leftatom,
|
---|
433 | current_bond->rightatom
|
---|
434 | };
|
---|
435 |
|
---|
436 | // remove the edge
|
---|
437 | #ifndef NDEBUG
|
---|
438 | const bool status =
|
---|
439 | #endif
|
---|
440 | BGcreator.removeEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
|
---|
441 | ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
|
---|
442 |
|
---|
443 | // gather nodes for either atom
|
---|
444 | BoostGraphHelpers::Nodeset_t bondside_set[MAX_sides];
|
---|
445 | BreadthFirstSearchGatherer::distance_map_t distance_map[MAX_sides];
|
---|
446 | for (size_t side=leftside;side<MAX_sides;++side) {
|
---|
447 | bondside_set[side] = NodeGatherer(bondatom[side]->getId(), max_distance);
|
---|
448 | distance_map[side] = NodeGatherer.getDistances();
|
---|
449 | std::sort(bondside_set[side].begin(), bondside_set[side].end());
|
---|
450 | }
|
---|
451 |
|
---|
452 | // re-add edge
|
---|
453 | BGcreator.addEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
|
---|
454 |
|
---|
455 | // do for both leftatom and rightatom of bond
|
---|
456 | for (size_t side = leftside; side < MAX_sides; ++side) {
|
---|
457 | const double &bondforce = projected_forces[0][side][index];
|
---|
458 | const double &oldbondforce = projected_forces[1][side][index];
|
---|
459 | const double bondforcedifference = fabs(bondforce - oldbondforce);
|
---|
460 | LOG(4, "DEBUG: bondforce for " << (side == leftside ? "left" : "right")
|
---|
461 | << " side of bond is " << bondforce);
|
---|
462 | LOG(4, "DEBUG: oldbondforce for " << (side == leftside ? "left" : "right")
|
---|
463 | << " side of bond is " << oldbondforce);
|
---|
464 | // if difference or bondforce itself is zero, do nothing
|
---|
465 | if ((fabs(bondforce) < MYEPSILON) || (fabs(bondforcedifference) < MYEPSILON))
|
---|
466 | continue;
|
---|
467 |
|
---|
468 | // get BondVector to bond
|
---|
469 | const BondVectors::mapped_t::const_iterator bviter =
|
---|
470 | bondvectors.find(current_bond);
|
---|
471 | ASSERT( bviter != bondvectors.end(),
|
---|
472 | "ForceAnnealing() - cannot find current_bond ?");
|
---|
473 | ASSERT( fabs(bviter->second.Norm() -1.) < MYEPSILON,
|
---|
474 | "ForceAnnealing() - norm of BondVector is not one");
|
---|
475 | const Vector &BondVector = bviter->second;
|
---|
476 |
|
---|
477 | // calculate gradient and position differences for stepwidth
|
---|
478 | const Vector currentGradient = bondforce * BondVector;
|
---|
479 | LOG(4, "DEBUG: current projected gradient for "
|
---|
480 | << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient);
|
---|
481 | const Vector &oldPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
|
---|
482 | const Vector ¤tPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0);
|
---|
483 | const Vector PositionDifference = currentPosition - oldPosition;
|
---|
484 | LOG(4, "DEBUG: old position is " << oldPosition);
|
---|
485 | LOG(4, "DEBUG: current position is " << currentPosition);
|
---|
486 | LOG(4, "DEBUG: difference in position is " << PositionDifference);
|
---|
487 | LOG(4, "DEBUG: bondvector is " << BondVector);
|
---|
488 | const double projected_PositionDifference = PositionDifference.ScalarProduct(BondVector);
|
---|
489 | LOG(4, "DEBUG: difference in position projected onto bondvector is "
|
---|
490 | << projected_PositionDifference);
|
---|
491 | LOG(4, "DEBUG: abs. difference in forces is " << bondforcedifference);
|
---|
492 |
|
---|
493 | // calculate step width
|
---|
494 | double stepwidth =
|
---|
495 | fabs(projected_PositionDifference)/bondforcedifference;
|
---|
496 | if (fabs(stepwidth) < 1e-10) {
|
---|
497 | // dont' warn in first step, deltat usage normal
|
---|
498 | if (currentStep != 1)
|
---|
499 | ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
|
---|
500 | stepwidth = currentDeltat;
|
---|
501 | }
|
---|
502 | Vector PositionUpdate = stepwidth * currentGradient;
|
---|
503 | LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
|
---|
504 |
|
---|
505 | // add PositionUpdate for all nodes in the bondside_set
|
---|
506 | for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set[side].begin();
|
---|
507 | setiter != bondside_set[side].end(); ++setiter) {
|
---|
508 | const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
|
---|
509 | = distance_map[side].find(*setiter);
|
---|
510 | ASSERT( diter != distance_map[side].end(),
|
---|
511 | "ForceAnnealing() - could not find distance to an atom.");
|
---|
512 | const double factor = pow(damping_factor, diter->second+1);
|
---|
513 | LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
|
---|
514 | << factor << "*" << PositionUpdate);
|
---|
515 | if (GatheredUpdates.count((*setiter))) {
|
---|
516 | GatheredUpdates[(*setiter)] += factor*PositionUpdate;
|
---|
517 | } else {
|
---|
518 | GatheredUpdates.insert(
|
---|
519 | std::make_pair(
|
---|
520 | (*setiter),
|
---|
521 | factor*PositionUpdate) );
|
---|
522 | }
|
---|
523 | }
|
---|
524 | }
|
---|
525 | }
|
---|
526 |
|
---|
527 | for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
|
---|
528 | iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
|
---|
529 | atom &walker = *(*iter);
|
---|
530 | // extract largest components for showing progress of annealing
|
---|
531 | const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep-1 : 0);
|
---|
532 | for(size_t i=0;i<NDIM;++i)
|
---|
533 | maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
|
---|
534 |
|
---|
535 | // reset force vector for next step except on final one
|
---|
536 | if (currentStep != maxSteps)
|
---|
537 | walker.setAtomicForce(zeroVec);
|
---|
538 | }
|
---|
539 |
|
---|
540 | // apply the gathered updates
|
---|
541 | for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
|
---|
542 | iter != GatheredUpdates.end(); ++iter) {
|
---|
543 | const atomId_t &atomid = iter->first;
|
---|
544 | const Vector &update = iter->second;
|
---|
545 | atom* const walker = World::getInstance().getAtom(AtomById(atomid));
|
---|
546 | ASSERT( walker != NULL,
|
---|
547 | "ForceAnnealing() - walker with id "+toString(atomid)+" has suddenly disappeared.");
|
---|
548 | LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
|
---|
549 | << ", namely " << *walker);
|
---|
550 | walker->setPosition(
|
---|
551 | walker->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0)
|
---|
552 | + update);
|
---|
553 | walker->setAtomicVelocity(update);
|
---|
554 | // walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
|
---|
555 | }
|
---|
556 | }
|
---|
557 |
|
---|
558 | /** Reset function to unset static entities and artificial velocities.
|
---|
559 | *
|
---|
560 | */
|
---|
561 | void reset()
|
---|
562 | {
|
---|
563 | currentDeltat = 0.;
|
---|
564 | currentStep = 0;
|
---|
565 | }
|
---|
566 |
|
---|
567 | private:
|
---|
568 | //!> contains the current step in relation to maxsteps
|
---|
569 | static size_t currentStep;
|
---|
570 | //!> contains the maximum number of steps, determines initial and final step with currentStep
|
---|
571 | size_t maxSteps;
|
---|
572 | static double currentDeltat;
|
---|
573 | //!> minimum deltat for internal while loop (adaptive step width)
|
---|
574 | static double MinimumDeltat;
|
---|
575 | //!> contains the maximum bond graph distance up to which shifts of a single atom are spread
|
---|
576 | const int max_distance;
|
---|
577 | //!> the shifted is dampened by this factor with the power of the bond graph distance to the shift causing atom
|
---|
578 | const double damping_factor;
|
---|
579 | };
|
---|
580 |
|
---|
581 | template <class T>
|
---|
582 | double ForceAnnealing<T>::currentDeltat = 0.;
|
---|
583 | template <class T>
|
---|
584 | size_t ForceAnnealing<T>::currentStep = 0;
|
---|
585 | template <class T>
|
---|
586 | double ForceAnnealing<T>::MinimumDeltat = 1e-8;
|
---|
587 |
|
---|
588 | #endif /* FORCEANNEALING_HPP_ */
|
---|