source: src/Dynamics/ForceAnnealing.hpp@ 1e49e6

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

Added reset() function to ForceAnnealing.

  • Property mode set to 100644
File size: 6.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 "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 * \param _atoms set of atoms to integrate
45 * \param _Deltat time step width in atomic units
46 * \param _IsAngstroem whether length units are in angstroem or bohr radii
47 * \param _maxSteps number of optimization steps to perform
48 */
49 ForceAnnealing(
50 AtomSetMixin<T> &_atoms,
51 double _Deltat,
52 bool _IsAngstroem,
53 const size_t _maxSteps) :
54 AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem),
55 maxSteps(_maxSteps)
56 {}
57 /** Destructor of class ForceAnnealing.
58 *
59 */
60 ~ForceAnnealing()
61 {}
62
63 /** Performs Gradient optimization.
64 *
65 * We assume that forces have just been calculated.
66 *
67 *
68 * \param NextStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
69 * \param offset offset in matrix file to the first force component
70 * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
71 */
72 void operator()(const int NextStep, const size_t offset)
73 {
74 // make sum of forces equal zero
75 AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(offset,NextStep);
76
77 // are we in initial step? Then set static entities
78 if (currentStep == 0) {
79 currentDeltat = AtomicForceManipulator<T>::Deltat;
80 currentStep = 1;
81 LOG(2, "DEBUG: Initial step, setting values, current step is #" << currentStep);
82 } else {
83 ++currentStep;
84 LOG(2, "DEBUG: current step is #" << currentStep);
85 }
86
87 Vector maxComponents(zeroVec);
88 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
89 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
90 // atom's force vector gives steepest descent direction
91 const Vector currentPosition = (*iter)->getPosition();
92 const Vector currentGradient = (*iter)->getAtomicForce();
93 LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
94
95 // artificial update: deltat may be considered as 1/2 s^2 units, mass
96 // is neglected deliberately as this makes all atoms equally fast or
97 // hydrogens slower (and they need to wait for other atoms to arrive at
98 // final position).
99 Vector PositionUpdate = currentDeltat * currentGradient;
100 LOG(3, "DEBUG: Update would be " << PositionUpdate);
101
102 // extract largest components for showing progress of annealing
103 for(size_t i=0;i<NDIM;++i)
104 if (currentGradient[i] > maxComponents[i])
105 maxComponents[i] = currentGradient[i];
106
107 // are we in initial step? Then don't check against velocity
108 if ((currentStep > 1) && (!(*iter)->getAtomicVelocity().IsZero()))
109 // update with currentDelta tells us how the current gradient relates to
110 // the last one: If it has become larger, reduce currentDelta
111 if ((PositionUpdate.ScalarProduct((*iter)->getAtomicVelocity()) < 0)
112 && (currentDeltat > MinimumDeltat)) {
113 currentDeltat = .5*currentDeltat;
114 LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate.NormSquared()
115 << " > " << (*iter)->getAtomicVelocity().NormSquared()
116 << ", decreasing deltat: " << currentDeltat);
117 PositionUpdate = currentDeltat * currentGradient;
118 }
119
120 // finally set new values
121 (*iter)->setPosition(currentPosition + PositionUpdate);
122 (*iter)->setAtomicVelocity(PositionUpdate);
123 //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
124// (*iter)->VelocityVerletUpdateU((*iter)->getId(), NextStep-1, Deltat, IsAngstroem);
125
126 // reset force vector for next step except on final one
127 if (currentStep != maxSteps)
128 (*iter)->setAtomicForce(zeroVec);
129 }
130
131 LOG(1, "STATUS: Largest remaining force components at step #"
132 << currentStep << " are " << maxComponents);
133
134 // are we in final step? Remember to reset static entities
135 if (currentStep == maxSteps) {
136 LOG(2, "DEBUG: Final step, resetting values");
137 reset();
138 }
139 }
140
141 /** Reset function to unset static entities and artificial velocities.
142 *
143 */
144 void reset()
145 {
146 currentDeltat = 0.;
147 currentStep = 0;
148
149 // reset (artifical) velocities
150 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
151 iter != AtomicForceManipulator<T>::atoms.end(); ++iter)
152 (*iter)->setAtomicVelocity(zeroVec);
153 }
154
155private:
156 //!> contains the current step in relation to maxsteps
157 static size_t currentStep;
158 //!> contains the maximum number of steps, determines initial and final step with currentStep
159 size_t maxSteps;
160 static double currentDeltat;
161 //!> minimum deltat for internal while loop (adaptive step width)
162 static double MinimumDeltat;
163};
164
165template <class T>
166double ForceAnnealing<T>::currentDeltat = 0.;
167template <class T>
168size_t ForceAnnealing<T>::currentStep = 0;
169template <class T>
170double ForceAnnealing<T>::MinimumDeltat = 1e-8;
171
172#endif /* FORCEANNEALING_HPP_ */
Note: See TracBrowser for help on using the repository browser.