source: src/Dynamics/ForceAnnealing.hpp@ 0c4f24

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

anneal() returns maxComponents, does not neet _offset.

  • Property mode set to 100644
File size: 24.0 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
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 */
56template <class T>
57class ForceAnnealing : public AtomicForceManipulator<T>
58{
59public:
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 maxComponents = anneal(_CurrentTimeStep);
117 } else {
118 ++currentStep;
119 LOG(2, "DEBUG: current step is #" << currentStep);
120
121 // bond graph annealing is always followed by a normal annealing
122 if (_UseBondgraph)
123 maxComponents = annealWithBondGraph(_CurrentTimeStep);
124 // cannot store RemnantGradient in Atom's Force as it ruins BB stepwidth calculation
125 else
126 maxComponents = anneal(_CurrentTimeStep);
127 }
128
129
130 LOG(1, "STATUS: Largest remaining force components at step #"
131 << currentStep << " are " << maxComponents);
132
133 // are we in final step? Remember to reset static entities
134 if (currentStep == maxSteps) {
135 LOG(2, "DEBUG: Final step, resetting values");
136 reset();
137 }
138 }
139
140 /** Helper function to calculate the Barzilai-Borwein stepwidth.
141 *
142 * \param _PositionDifference difference in position between current and last step
143 * \param _GradientDifference difference in gradient between current and last step
144 * \return step width according to Barzilai-Borwein
145 */
146 double getBarzilaiBorweinStepwidth(const Vector &_PositionDifference, const Vector &_GradientDifference)
147 {
148 double stepwidth = 0.;
149 if (_GradientDifference.NormSquared() > MYEPSILON)
150 stepwidth = fabs(_PositionDifference.ScalarProduct(_GradientDifference))/
151 _GradientDifference.NormSquared();
152 if (fabs(stepwidth) < 1e-10) {
153 // dont' warn in first step, deltat usage normal
154 if (currentStep != 1)
155 ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
156 stepwidth = currentDeltat;
157 }
158 return stepwidth;
159 }
160
161 /** Performs Gradient optimization on the atoms.
162 *
163 * We assume that forces have just been calculated.
164 *
165 * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
166 * \return to be filled with maximum force component over all atoms
167 */
168 Vector anneal(
169 const int CurrentTimeStep)
170 {
171 Vector maxComponents;
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 return maxComponents;
219 }
220
221 // knowing the number of bonds in total, we can setup the storage for the
222 // projected forces
223 enum whichatom_t {
224 leftside=0,
225 rightside=1,
226 MAX_sides
227 };
228
229 /** Helper function to put bond force into a container.
230 *
231 * \param _walker atom
232 * \param _current_bond current bond of \a _walker
233 * \param _timestep time step
234 * \param _force calculated bond force
235 * \param _bv bondvectors for obtaining the correct index
236 * \param _projected_forces container
237 */
238 static void ForceStore(
239 const atom &_walker,
240 const bond::ptr &_current_bond,
241 const size_t &_timestep,
242 const double _force,
243 const BondVectors &_bv,
244 std::vector< // time step
245 std::vector< // which bond side
246 std::vector<double> > // over all bonds
247 > &_projected_forces)
248 {
249 std::vector<double> &forcelist = (&_walker == _current_bond->leftatom) ?
250 _projected_forces[_timestep][leftside] : _projected_forces[_timestep][rightside];
251 const size_t index = _bv.getIndexForBond(_current_bond);
252 ASSERT( index != (size_t)-1,
253 "ForceAnnealing() - could not find bond "+toString(*_current_bond)
254 +" in bondvectors");
255 forcelist[index] = _force;
256 }
257
258 /** Performs Gradient optimization on the bonds.
259 *
260 * We assume that forces have just been calculated. These forces are projected
261 * onto the bonds and these are annealed subsequently by moving atoms in the
262 * bond neighborhood on either side conjunctively.
263 *
264 *
265 * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
266 * \param maxComponents to be filled with maximum force component over all atoms
267 */
268 Vector annealWithBondGraph(
269 const int CurrentTimeStep)
270 {
271 Vector maxComponents;
272
273 // get nodes on either side of selected bond via BFS discovery
274 BoostGraphCreator BGcreator;
275 BGcreator.createFromRange(
276 AtomicForceManipulator<T>::atoms.begin(),
277 AtomicForceManipulator<T>::atoms.end(),
278 AtomicForceManipulator<T>::atoms.size(),
279 BreadthFirstSearchGatherer::AlwaysTruePredicate);
280 BreadthFirstSearchGatherer NodeGatherer(BGcreator);
281
282 /// We assume that a force is local, i.e. a bond is too short yet and hence
283 /// the atom needs to be moved. However, all the adjacent (bound) atoms might
284 /// already be at the perfect distance. If we just move the atom alone, we ruin
285 /// all the other bonds. Hence, it would be sensible to move every atom found
286 /// through the bond graph in the direction of the force as well by the same
287 /// PositionUpdate. This is almost what we are going to do.
288
289 /// One issue is: If we need to shorten bond, then we use the PositionUpdate
290 /// also on the the other bond partner already. This is because it is in the
291 /// direction of the bond. Therefore, the update is actually performed twice on
292 /// each bond partner, i.e. the step size is twice as large as it should be.
293 /// This problem only occurs when bonds need to be shortened, not when they
294 /// need to be made longer (then the force vector is facing the other
295 /// direction than the bond vector).
296 /// As a remedy we need to average the force on either end of the bond and
297 /// check whether each gradient points inwards out or outwards with respect
298 /// to the bond and then shift accordingly.
299
300 /// One more issue is that the projection onto the bond directions does not
301 /// recover the gradient but may be larger as the bond directions are a
302 /// generating system and not a basis (e.g. 3 bonds on a plane where 2 would
303 /// suffice to span the plane). To this end, we need to account for the
304 /// overestimation and obtain a weighting for each bond.
305
306 // initialize helper class for bond vectors using bonds from range of atoms
307 BondVectors bv;
308 bv.setFromAtomRange< T >(
309 AtomicForceManipulator<T>::atoms.begin(),
310 AtomicForceManipulator<T>::atoms.end(),
311 CurrentTimeStep);
312 const BondVectors::container_t &sorted_bonds = bv.getSorted();
313
314 std::vector< // time step
315 std::vector< // which bond side
316 std::vector<double> > // over all bonds
317 > projected_forces(2); // one for leftatoms, one for rightatoms (and for both time steps)
318 for (size_t i=0;i<2;++i) {
319 projected_forces[i].resize(MAX_sides);
320 for (size_t j=0;j<MAX_sides;++j)
321 projected_forces[i][j].resize(sorted_bonds.size(), 0.);
322 }
323
324 // for each atom we need to gather weights and then project the gradient
325 typedef std::map<atomId_t, BondVectors::weights_t > weights_per_atom_t;
326 std::vector<weights_per_atom_t> weights_per_atom(2);
327 typedef std::map<atomId_t, Vector> RemnantGradient_per_atom_t;
328 RemnantGradient_per_atom_t RemnantGradient_per_atom;
329 for (size_t timestep = 0; timestep <= 1; ++timestep) {
330 const size_t CurrentStep = CurrentTimeStep-timestep-1 >= 0 ? CurrentTimeStep-timestep-1 : 0;
331 LOG(2, "DEBUG: CurrentTimeStep is " << CurrentTimeStep
332 << ", timestep is " << timestep
333 << ", and CurrentStep is " << CurrentStep);
334
335 for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin();
336 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
337 const atom &walker = *(*iter);
338 const Vector &walkerGradient = walker.getAtomicForceAtStep(CurrentStep);
339 LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
340 << walker << " is " << walkerGradient << " with magnitude of "
341 << walkerGradient.Norm());
342
343 const BondList& ListOfBonds = walker.getListOfBonds();
344 if (walkerGradient.Norm() > MYEPSILON) {
345
346 // gather subset of BondVectors for the current atom
347 const std::vector<Vector> BondVectors =
348 bv.getAtomsBondVectorsAtStep(walker, CurrentStep);
349
350 // go through all its bonds and calculate what magnitude is represented
351 // by the others i.e. sum of scalar products against other bonds
352 const std::pair<weights_per_atom_t::iterator, bool> inserter =
353 weights_per_atom[timestep].insert(
354 std::make_pair(walker.getId(),
355 bv.getWeightsForAtomAtStep(walker, BondVectors, CurrentStep)) );
356 ASSERT( inserter.second,
357 "ForceAnnealing::operator() - weight map for atom "+toString(walker)
358 +" and time step "+toString(timestep)+" already filled?");
359 BondVectors::weights_t &weights = inserter.first->second;
360 ASSERT( weights.size() == ListOfBonds.size(),
361 "ForceAnnealing::operator() - number of weights "
362 +toString(weights.size())+" does not match number of bonds "
363 +toString(ListOfBonds.size())+", error in calculation?");
364
365 // projected gradient over all bonds and place in one of projected_forces
366 // using the obtained weights
367 BondVectors::forcestore_t forcestoring =
368 boost::bind(&ForceAnnealing::ForceStore, _1, _2, _3, _4,
369 boost::cref(bv), boost::ref(projected_forces));
370 const Vector RemnantGradient = bv.getRemnantGradientForAtomAtStep(
371 walker, walkerGradient, BondVectors, weights, timestep, forcestoring
372 );
373 RemnantGradient_per_atom.insert( std::make_pair(walker.getId(), RemnantGradient) );
374 } else {
375 LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
376 << MYEPSILON << " for atom " << walker);
377 // note that projected_forces is initialized to full length and filled
378 // with zeros. Hence, nothing to do here
379 }
380 }
381 }
382
383 // step through each bond and shift the atoms
384 std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
385
386 LOG(3, "DEBUG: current step is " << currentStep << ", given time step is " << CurrentTimeStep);
387 const BondVectors::mapped_t bondvectors = bv.getBondVectorsAtStep(CurrentTimeStep);
388
389 for (BondVectors::container_t::const_iterator bondsiter = sorted_bonds.begin();
390 bondsiter != sorted_bonds.end(); ++bondsiter) {
391 const bond::ptr &current_bond = *bondsiter;
392 const size_t index = bv.getIndexForBond(current_bond);
393 const atom* bondatom[MAX_sides] = {
394 current_bond->leftatom,
395 current_bond->rightatom
396 };
397
398 // remove the edge
399#ifndef NDEBUG
400 const bool status =
401#endif
402 BGcreator.removeEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
403 ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
404
405 // gather nodes for either atom
406 BoostGraphHelpers::Nodeset_t bondside_set[MAX_sides];
407 BreadthFirstSearchGatherer::distance_map_t distance_map[MAX_sides];
408 for (size_t side=leftside;side<MAX_sides;++side) {
409 bondside_set[side] = NodeGatherer(bondatom[side]->getId(), max_distance);
410 distance_map[side] = NodeGatherer.getDistances();
411 std::sort(bondside_set[side].begin(), bondside_set[side].end());
412 }
413
414 // re-add edge
415 BGcreator.addEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
416
417 // do for both leftatom and rightatom of bond
418 for (size_t side = leftside; side < MAX_sides; ++side) {
419 const double &bondforce = projected_forces[0][side][index];
420 const double &oldbondforce = projected_forces[1][side][index];
421 const double bondforcedifference = fabs(bondforce - oldbondforce);
422 LOG(4, "DEBUG: bondforce for " << (side == leftside ? "left" : "right")
423 << " side of bond is " << bondforce);
424 LOG(4, "DEBUG: oldbondforce for " << (side == leftside ? "left" : "right")
425 << " side of bond is " << oldbondforce);
426 // if difference or bondforce itself is zero, do nothing
427 if ((fabs(bondforce) < MYEPSILON) || (fabs(bondforcedifference) < MYEPSILON))
428 continue;
429
430 // get BondVector to bond
431 const BondVectors::mapped_t::const_iterator bviter =
432 bondvectors.find(current_bond);
433 ASSERT( bviter != bondvectors.end(),
434 "ForceAnnealing() - cannot find current_bond ?");
435 ASSERT( fabs(bviter->second.Norm() -1.) < MYEPSILON,
436 "ForceAnnealing() - norm of BondVector is not one");
437 const Vector &BondVector = bviter->second;
438
439 // calculate gradient and position differences for stepwidth
440 const Vector currentGradient = bondforce * BondVector;
441 LOG(4, "DEBUG: current projected gradient for "
442 << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient);
443 const Vector &oldPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
444 const Vector &currentPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0);
445 const Vector PositionDifference = currentPosition - oldPosition;
446 LOG(4, "DEBUG: old position is " << oldPosition);
447 LOG(4, "DEBUG: current position is " << currentPosition);
448 LOG(4, "DEBUG: difference in position is " << PositionDifference);
449 LOG(4, "DEBUG: bondvector is " << BondVector);
450 const double projected_PositionDifference = PositionDifference.ScalarProduct(BondVector);
451 LOG(4, "DEBUG: difference in position projected onto bondvector is "
452 << projected_PositionDifference);
453 LOG(4, "DEBUG: abs. difference in forces is " << bondforcedifference);
454
455 // calculate step width
456 double stepwidth =
457 fabs(projected_PositionDifference)/bondforcedifference;
458 if (fabs(stepwidth) < 1e-10) {
459 // dont' warn in first step, deltat usage normal
460 if (currentStep != 1)
461 ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
462 stepwidth = currentDeltat;
463 }
464 Vector PositionUpdate = stepwidth * currentGradient;
465 LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
466
467 // add PositionUpdate for all nodes in the bondside_set
468 for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set[side].begin();
469 setiter != bondside_set[side].end(); ++setiter) {
470 const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
471 = distance_map[side].find(*setiter);
472 ASSERT( diter != distance_map[side].end(),
473 "ForceAnnealing() - could not find distance to an atom.");
474 const double factor = pow(damping_factor, diter->second+1);
475 LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
476 << factor << "*" << PositionUpdate);
477 if (GatheredUpdates.count((*setiter))) {
478 GatheredUpdates[(*setiter)] += factor*PositionUpdate;
479 } else {
480 GatheredUpdates.insert(
481 std::make_pair(
482 (*setiter),
483 factor*PositionUpdate) );
484 }
485 }
486 }
487 }
488
489 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
490 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
491 atom &walker = *(*iter);
492 // extract largest components for showing progress of annealing
493 const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep-1 : 0);
494 for(size_t i=0;i<NDIM;++i)
495 maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
496
497 // reset force vector for next step except on final one
498 if (currentStep != maxSteps)
499 walker.setAtomicForce(zeroVec);
500 }
501
502 // remove center of weight translation from gathered updates
503 Vector CommonTranslation;
504 for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
505 iter != GatheredUpdates.end(); ++iter) {
506 const Vector &update = iter->second;
507 CommonTranslation += update;
508 }
509 CommonTranslation *= 1./(double)GatheredUpdates.size();
510 LOG(3, "DEBUG: Subtracting common translation " << CommonTranslation
511 << " from all updates.");
512
513 // apply the gathered updates and set remnant gradients for atomic annealing
514 for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
515 iter != GatheredUpdates.end(); ++iter) {
516 const atomId_t &atomid = iter->first;
517 const Vector &update = iter->second;
518 atom* const walker = World::getInstance().getAtom(AtomById(atomid));
519 ASSERT( walker != NULL,
520 "ForceAnnealing() - walker with id "+toString(atomid)+" has suddenly disappeared.");
521 LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
522 << ", namely " << *walker);
523 walker->setPosition(
524 walker->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0)
525 + update - CommonTranslation);
526 walker->setAtomicVelocity(update);
527// walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
528 }
529
530 return maxComponents;
531 }
532
533 /** Reset function to unset static entities and artificial velocities.
534 *
535 */
536 void reset()
537 {
538 currentDeltat = 0.;
539 currentStep = 0;
540 }
541
542private:
543 //!> contains the current step in relation to maxsteps
544 static size_t currentStep;
545 //!> contains the maximum number of steps, determines initial and final step with currentStep
546 size_t maxSteps;
547 static double currentDeltat;
548 //!> minimum deltat for internal while loop (adaptive step width)
549 static double MinimumDeltat;
550 //!> contains the maximum bond graph distance up to which shifts of a single atom are spread
551 const int max_distance;
552 //!> the shifted is dampened by this factor with the power of the bond graph distance to the shift causing atom
553 const double damping_factor;
554};
555
556template <class T>
557double ForceAnnealing<T>::currentDeltat = 0.;
558template <class T>
559size_t ForceAnnealing<T>::currentStep = 0;
560template <class T>
561double ForceAnnealing<T>::MinimumDeltat = 1e-8;
562
563#endif /* FORCEANNEALING_HPP_ */
Note: See TracBrowser for help on using the repository browser.