source: src/Dynamics/ForceAnnealing.hpp@ 825d33

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 825d33 was 825d33, checked in by Frederik Heber <frederik.heber@…>, 7 years ago

Extracted calculation of weights per atom into BondVectors.

  • Property mode set to 100644
File size: 23.4 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 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::map<atomId_t, BondVectors::weights_t > weights_per_atom_t;
295 std::vector<weights_per_atom_t> weights_per_atom(2);
296 for (size_t timestep = 0; timestep <= 1; ++timestep) {
297 const size_t CurrentStep = CurrentTimeStep-timestep-1 >= 0 ? CurrentTimeStep-timestep-1 : 0;
298 LOG(2, "DEBUG: CurrentTimeStep is " << CurrentTimeStep
299 << ", timestep is " << timestep
300 << ", and CurrentStep is " << CurrentStep);
301
302 for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin();
303 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
304 const atom &walker = *(*iter);
305 const Vector &walkerGradient = walker.getAtomicForceAtStep(CurrentStep);
306 LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
307 << walker << " is " << walkerGradient << " with magnitude of "
308 << walkerGradient.Norm());
309
310 const BondList& ListOfBonds = walker.getListOfBonds();
311 if (walkerGradient.Norm() > MYEPSILON) {
312
313 // gather subset of BondVectors for the current atom
314 const std::vector<Vector> BondVectors =
315 bv.getAtomsBondVectorsAtStep(walker, CurrentStep);
316
317 // go through all its bonds and calculate what magnitude is represented
318 // by the others i.e. sum of scalar products against other bonds
319 const std::pair<weights_per_atom_t::iterator, bool> inserter =
320 weights_per_atom[timestep].insert(
321 std::make_pair(walker.getId(),
322 bv.getWeightsForAtomAtStep(walker, CurrentStep)) );
323 ASSERT( inserter.second,
324 "ForceAnnealing::operator() - weight map for atom "+toString(walker)
325 +" and time step "+toString(timestep)+" already filled?");
326 BondVectors::weights_t &weights = inserter.first->second;
327 ASSERT( weights.size() == ListOfBonds.size(),
328 "ForceAnnealing::operator() - number of weights "
329 +toString(weights.size())+" does not match number of bonds "
330 +toString(ListOfBonds.size())+", error in calculation?");
331
332 // projected gradient over all bonds and place in one of projected_forces
333 // using the obtained weights
334 {
335 BondVectors::weights_t::const_iterator weightiter = weights.begin();
336 std::vector<Vector>::const_iterator vectoriter = BondVectors.begin();
337 Vector forcesum_check;
338 for(BondList::const_iterator bonditer = ListOfBonds.begin();
339 bonditer != ListOfBonds.end(); ++bonditer, ++weightiter, ++vectoriter) {
340 const bond::ptr &current_bond = *bonditer;
341 const Vector &BondVector = *vectoriter;
342
343 std::vector<double> &forcelist = (&walker == current_bond->leftatom) ?
344 projected_forces[timestep][leftside] : projected_forces[timestep][rightside];
345 const size_t index = bv.getIndexForBond(current_bond);
346 ASSERT( index != (size_t)-1,
347 "ForceAnnealing() - could not find bond "+toString(*current_bond)
348 +" in bondvectors");
349 forcelist[index] = (*weightiter)*walkerGradient.ScalarProduct(BondVector);
350 LOG(4, "DEBUG: BondVector " << BondVector << " receives projected force of "
351 << forcelist[index]);
352 forcesum_check += forcelist[index] * BondVector;
353 }
354 ASSERT( weightiter == weights.end(),
355 "ForceAnnealing - weightiter is not at end when it should be.");
356 ASSERT( vectoriter == BondVectors.end(),
357 "ForceAnnealing - vectoriter is not at end when it should be.");
358 LOG(3, "DEBUG: sum of projected forces is " << forcesum_check);
359 }
360
361 } else {
362 LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
363 << MYEPSILON << " for atom " << walker);
364 // note that projected_forces is initialized to full length and filled
365 // with zeros. Hence, nothing to do here
366 }
367 }
368 }
369
370 // step through each bond and shift the atoms
371 std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
372
373 LOG(3, "DEBUG: current step is " << currentStep << ", given time step is " << CurrentTimeStep);
374 const BondVectors::mapped_t bondvectors = bv.getBondVectorsAtStep(CurrentTimeStep);
375
376 for (BondVectors::container_t::const_iterator bondsiter = sorted_bonds.begin();
377 bondsiter != sorted_bonds.end(); ++bondsiter) {
378 const bond::ptr &current_bond = *bondsiter;
379 const size_t index = bv.getIndexForBond(current_bond);
380 const atom* bondatom[MAX_sides] = {
381 current_bond->leftatom,
382 current_bond->rightatom
383 };
384
385 // remove the edge
386#ifndef NDEBUG
387 const bool status =
388#endif
389 BGcreator.removeEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
390 ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
391
392 // gather nodes for either atom
393 BoostGraphHelpers::Nodeset_t bondside_set[MAX_sides];
394 BreadthFirstSearchGatherer::distance_map_t distance_map[MAX_sides];
395 for (size_t side=leftside;side<MAX_sides;++side) {
396 bondside_set[side] = NodeGatherer(bondatom[side]->getId(), max_distance);
397 distance_map[side] = NodeGatherer.getDistances();
398 std::sort(bondside_set[side].begin(), bondside_set[side].end());
399 }
400
401 // re-add edge
402 BGcreator.addEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
403
404 // do for both leftatom and rightatom of bond
405 for (size_t side = leftside; side < MAX_sides; ++side) {
406 const double &bondforce = projected_forces[0][side][index];
407 const double &oldbondforce = projected_forces[1][side][index];
408 const double bondforcedifference = fabs(bondforce - oldbondforce);
409 LOG(4, "DEBUG: bondforce for " << (side == leftside ? "left" : "right")
410 << " side of bond is " << bondforce);
411 LOG(4, "DEBUG: oldbondforce for " << (side == leftside ? "left" : "right")
412 << " side of bond is " << oldbondforce);
413 // if difference or bondforce itself is zero, do nothing
414 if ((fabs(bondforce) < MYEPSILON) || (fabs(bondforcedifference) < MYEPSILON))
415 continue;
416
417 // get BondVector to bond
418 const BondVectors::mapped_t::const_iterator bviter =
419 bondvectors.find(current_bond);
420 ASSERT( bviter != bondvectors.end(),
421 "ForceAnnealing() - cannot find current_bond ?");
422 ASSERT( fabs(bviter->second.Norm() -1.) < MYEPSILON,
423 "ForceAnnealing() - norm of BondVector is not one");
424 const Vector &BondVector = bviter->second;
425
426 // calculate gradient and position differences for stepwidth
427 const Vector currentGradient = bondforce * BondVector;
428 LOG(4, "DEBUG: current projected gradient for "
429 << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient);
430 const Vector &oldPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
431 const Vector &currentPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0);
432 const Vector PositionDifference = currentPosition - oldPosition;
433 LOG(4, "DEBUG: old position is " << oldPosition);
434 LOG(4, "DEBUG: current position is " << currentPosition);
435 LOG(4, "DEBUG: difference in position is " << PositionDifference);
436 LOG(4, "DEBUG: bondvector is " << BondVector);
437 const double projected_PositionDifference = PositionDifference.ScalarProduct(BondVector);
438 LOG(4, "DEBUG: difference in position projected onto bondvector is "
439 << projected_PositionDifference);
440 LOG(4, "DEBUG: abs. difference in forces is " << bondforcedifference);
441
442 // calculate step width
443 double stepwidth =
444 fabs(projected_PositionDifference)/bondforcedifference;
445 if (fabs(stepwidth) < 1e-10) {
446 // dont' warn in first step, deltat usage normal
447 if (currentStep != 1)
448 ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
449 stepwidth = currentDeltat;
450 }
451 Vector PositionUpdate = stepwidth * currentGradient;
452 LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
453
454 // add PositionUpdate for all nodes in the bondside_set
455 for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set[side].begin();
456 setiter != bondside_set[side].end(); ++setiter) {
457 const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
458 = distance_map[side].find(*setiter);
459 ASSERT( diter != distance_map[side].end(),
460 "ForceAnnealing() - could not find distance to an atom.");
461 const double factor = pow(damping_factor, diter->second+1);
462 LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
463 << factor << "*" << PositionUpdate);
464 if (GatheredUpdates.count((*setiter))) {
465 GatheredUpdates[(*setiter)] += factor*PositionUpdate;
466 } else {
467 GatheredUpdates.insert(
468 std::make_pair(
469 (*setiter),
470 factor*PositionUpdate) );
471 }
472 }
473 }
474 }
475
476 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
477 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
478 atom &walker = *(*iter);
479 // extract largest components for showing progress of annealing
480 const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep-1 : 0);
481 for(size_t i=0;i<NDIM;++i)
482 maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
483
484 // reset force vector for next step except on final one
485 if (currentStep != maxSteps)
486 walker.setAtomicForce(zeroVec);
487 }
488
489 // apply the gathered updates
490 for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
491 iter != GatheredUpdates.end(); ++iter) {
492 const atomId_t &atomid = iter->first;
493 const Vector &update = iter->second;
494 atom* const walker = World::getInstance().getAtom(AtomById(atomid));
495 ASSERT( walker != NULL,
496 "ForceAnnealing() - walker with id "+toString(atomid)+" has suddenly disappeared.");
497 LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
498 << ", namely " << *walker);
499 walker->setPosition(
500 walker->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0)
501 + update);
502 walker->setAtomicVelocity(update);
503// walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
504 }
505 }
506
507 /** Reset function to unset static entities and artificial velocities.
508 *
509 */
510 void reset()
511 {
512 currentDeltat = 0.;
513 currentStep = 0;
514 }
515
516private:
517 //!> contains the current step in relation to maxsteps
518 static size_t currentStep;
519 //!> contains the maximum number of steps, determines initial and final step with currentStep
520 size_t maxSteps;
521 static double currentDeltat;
522 //!> minimum deltat for internal while loop (adaptive step width)
523 static double MinimumDeltat;
524 //!> contains the maximum bond graph distance up to which shifts of a single atom are spread
525 const int max_distance;
526 //!> the shifted is dampened by this factor with the power of the bond graph distance to the shift causing atom
527 const double damping_factor;
528};
529
530template <class T>
531double ForceAnnealing<T>::currentDeltat = 0.;
532template <class T>
533size_t ForceAnnealing<T>::currentStep = 0;
534template <class T>
535double ForceAnnealing<T>::MinimumDeltat = 1e-8;
536
537#endif /* FORCEANNEALING_HPP_ */
Note: See TracBrowser for help on using the repository browser.