source: src/Dynamics/ForceAnnealing.hpp@ 9bb8c8

AutomationFragmentation_failures Candidate_v1.6.1 ChemicalSpaceEvaluator Enhanced_StructuralOptimization_continued 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 9bb8c8 was 9bb8c8, checked in by Frederik Heber <frederik.heber@…>, 7 years ago

FIX: maxComponents was not chosen by absolute magnitude.

  • Property mode set to 100644
File size: 15.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 "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 "Descriptors/AtomIdDescriptor.hpp"
23#include "Dynamics/AtomicForceManipulator.hpp"
24#include "Fragmentation/ForceMatrix.hpp"
25#include "Graph/BoostGraphCreator.hpp"
26#include "Graph/BoostGraphHelpers.hpp"
27#include "Graph/BreadthFirstSearchGatherer.hpp"
28#include "Helpers/helpers.hpp"
29#include "Helpers/defs.hpp"
30#include "LinearAlgebra/Vector.hpp"
31#include "Thermostats/ThermoStatContainer.hpp"
32#include "Thermostats/Thermostat.hpp"
33#include "World.hpp"
34
35/** This class is the essential build block for performing structural optimization.
36 *
37 * Sadly, we have to use some static instances as so far values cannot be passed
38 * between actions. Hence, we need to store the current step and the adaptive-
39 * step width (we cannot perform a line search, as we have no control over the
40 * calculation of the forces).
41 *
42 * However, we do use the bond graph, i.e. if a single atom needs to be shifted
43 * to the left, then the whole molecule left of it is shifted, too. This is
44 * controlled by the \a max_distance parameter.
45 */
46template <class T>
47class ForceAnnealing : public AtomicForceManipulator<T>
48{
49public:
50 /** Constructor of class ForceAnnealing.
51 *
52 * \note We use a fixed delta t of 1.
53 *
54 * \param _atoms set of atoms to integrate
55 * \param _Deltat time step width in atomic units
56 * \param _IsAngstroem whether length units are in angstroem or bohr radii
57 * \param _maxSteps number of optimization steps to perform
58 * \param _max_distance up to this bond order is bond graph taken into account.
59 */
60 ForceAnnealing(
61 AtomSetMixin<T> &_atoms,
62 const double _Deltat,
63 bool _IsAngstroem,
64 const size_t _maxSteps,
65 const int _max_distance,
66 const double _damping_factor) :
67 AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem),
68 maxSteps(_maxSteps),
69 max_distance(_max_distance),
70 damping_factor(_damping_factor)
71 {}
72
73 /** Destructor of class ForceAnnealing.
74 *
75 */
76 ~ForceAnnealing()
77 {}
78
79 /** Performs Gradient optimization.
80 *
81 * We assume that forces have just been calculated.
82 *
83 *
84 * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
85 * \param offset offset in matrix file to the first force component
86 * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
87 */
88 void operator()(
89 const int _CurrentTimeStep,
90 const size_t _offset,
91 const bool _UseBondgraph)
92 {
93 // make sum of forces equal zero
94 AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(_offset, _CurrentTimeStep);
95
96 // are we in initial step? Then set static entities
97 Vector maxComponents(zeroVec);
98 if (currentStep == 0) {
99 currentDeltat = AtomicForceManipulator<T>::Deltat;
100 currentStep = 1;
101 LOG(2, "DEBUG: Initial step, setting values, current step is #" << currentStep);
102
103 // always use atomic annealing on first step
104 anneal(_CurrentTimeStep, _offset, maxComponents);
105 } else {
106 ++currentStep;
107 LOG(2, "DEBUG: current step is #" << currentStep);
108
109 if (_UseBondgraph)
110 annealWithBondGraph(_CurrentTimeStep, _offset, maxComponents);
111 else
112 anneal(_CurrentTimeStep, _offset, maxComponents);
113 }
114
115 LOG(1, "STATUS: Largest remaining force components at step #"
116 << currentStep << " are " << maxComponents);
117
118 // are we in final step? Remember to reset static entities
119 if (currentStep == maxSteps) {
120 LOG(2, "DEBUG: Final step, resetting values");
121 reset();
122 }
123 }
124
125 /** Performs Gradient optimization on the atoms.
126 *
127 * We assume that forces have just been calculated.
128 *
129 * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
130 * \param offset offset in matrix file to the first force component
131 * \param maxComponents to be filled with maximum force component over all atoms
132 */
133 void anneal(
134 const int CurrentTimeStep,
135 const size_t offset,
136 Vector &maxComponents)
137 {
138 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
139 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
140 // atom's force vector gives steepest descent direction
141 const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
142 const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
143 const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
144 const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
145 LOG(4, "DEBUG: oldPosition for atom " << **iter << " is " << oldPosition);
146 LOG(4, "DEBUG: currentPosition for atom " << **iter << " is " << currentPosition);
147 LOG(4, "DEBUG: oldGradient for atom " << **iter << " is " << oldGradient);
148 LOG(4, "DEBUG: currentGradient for atom " << **iter << " is " << currentGradient);
149// LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
150
151 // we use Barzilai-Borwein update with position reversed to get descent
152 const Vector PositionDifference = currentPosition - oldPosition;
153 const Vector GradientDifference = (currentGradient - oldGradient);
154 double stepwidth = 0.;
155 if (GradientDifference.Norm() > MYEPSILON)
156 stepwidth = fabs(PositionDifference.ScalarProduct(GradientDifference))/
157 GradientDifference.NormSquared();
158 if (fabs(stepwidth) < 1e-10) {
159 // dont' warn in first step, deltat usage normal
160 if (currentStep != 1)
161 ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
162 stepwidth = currentDeltat;
163 }
164 Vector PositionUpdate = stepwidth * currentGradient;
165 LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
166
167 // extract largest components for showing progress of annealing
168 for(size_t i=0;i<NDIM;++i)
169 maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
170
171 // are we in initial step? Then don't check against velocity
172 if ((currentStep > 1) && (!(*iter)->getAtomicVelocity().IsZero()))
173 // update with currentDelta tells us how the current gradient relates to
174 // the last one: If it has become larger, reduce currentDelta
175 if ((PositionUpdate.ScalarProduct((*iter)->getAtomicVelocity()) < 0)
176 && (currentDeltat > MinimumDeltat)) {
177 currentDeltat = .5*currentDeltat;
178 LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate.NormSquared()
179 << " > " << (*iter)->getAtomicVelocity().NormSquared()
180 << ", decreasing deltat: " << currentDeltat);
181 PositionUpdate = currentDeltat * currentGradient;
182 }
183 // finally set new values
184 (*iter)->setPosition(currentPosition + PositionUpdate);
185 (*iter)->setAtomicVelocity(PositionUpdate);
186 //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
187// (*iter)->VelocityVerletUpdateU((*iter)->getId(), CurrentTimeStep-1, Deltat, IsAngstroem);
188 }
189 }
190
191 /** Performs Gradient optimization on the bonds.
192 *
193 * We assume that forces have just been calculated. These forces are projected
194 * onto the bonds and these are annealed subsequently by moving atoms in the
195 * bond neighborhood on either side conjunctively.
196 *
197 *
198 * \param CurrentTimeStep current time step (i.e. t where \f$ t + \Delta t \f$ is in the sense of the velocity verlet)
199 * \param offset offset in matrix file to the first force component
200 * \param maxComponents to be filled with maximum force component over all atoms
201 */
202 void annealWithBondGraph(
203 const int CurrentTimeStep,
204 const size_t offset,
205 Vector &maxComponents)
206 {
207 // get nodes on either side of selected bond via BFS discovery
208// std::vector<atomId_t> atomids;
209// for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
210// iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
211// atomids.push_back((*iter)->getId());
212// }
213// ASSERT( atomids.size() == AtomicForceManipulator<T>::atoms.size(),
214// "ForceAnnealing() - could not gather all atomic ids?");
215 BoostGraphCreator BGcreator;
216 BGcreator.createFromRange(
217 AtomicForceManipulator<T>::atoms.begin(),
218 AtomicForceManipulator<T>::atoms.end(),
219 AtomicForceManipulator<T>::atoms.size(),
220 BreadthFirstSearchGatherer::AlwaysTruePredicate);
221 BreadthFirstSearchGatherer NodeGatherer(BGcreator);
222
223 std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
224 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
225 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
226 // atom's force vector gives steepest descent direction
227 const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
228 const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
229 const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
230 const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
231 LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
232
233 // we use Barzilai-Borwein update with position reversed to get descent
234 const Vector GradientDifference = (currentGradient - oldGradient);
235 const double stepwidth =
236 fabs((currentPosition - oldPosition).ScalarProduct(GradientDifference))/
237 GradientDifference.NormSquared();
238 Vector PositionUpdate = stepwidth * currentGradient;
239 if (fabs(stepwidth) < 1e-10) {
240 // dont' warn in first step, deltat usage normal
241 if (currentStep != 1)
242 ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
243 PositionUpdate = currentDeltat * currentGradient;
244 }
245 LOG(3, "DEBUG: Update would be " << PositionUpdate);
246
247// // add update to central atom
248// const atomId_t atomid = (*iter)->getId();
249// if (GatheredUpdates.count(atomid)) {
250// GatheredUpdates[atomid] += PositionUpdate;
251// } else
252// GatheredUpdates.insert( std::make_pair(atomid, PositionUpdate) );
253
254 // We assume that a force is local, i.e. a bond is too short yet and hence
255 // the atom needs to be moved. However, all the adjacent (bound) atoms might
256 // already be at the perfect distance. If we just move the atom alone, we ruin
257 // all the other bonds. Hence, it would be sensible to move every atom found
258 // through the bond graph in the direction of the force as well by the same
259 // PositionUpdate. This is just what we are going to do.
260
261 /// get all nodes from bonds in the direction of the current force
262
263 // remove edges facing in the wrong direction
264 std::vector<bond::ptr> removed_bonds;
265 const BondList& ListOfBonds = (*iter)->getListOfBonds();
266 for(BondList::const_iterator bonditer = ListOfBonds.begin();
267 bonditer != ListOfBonds.end(); ++bonditer) {
268 const bond &current_bond = *(*bonditer);
269 LOG(2, "DEBUG: Looking at bond " << current_bond);
270 Vector BondVector = (*iter)->getPositionAtStep(CurrentTimeStep);
271 BondVector -= ((*iter)->getId() == current_bond.rightatom->getId())
272 ? current_bond.rightatom->getPositionAtStep(CurrentTimeStep) : current_bond.leftatom->getPositionAtStep(CurrentTimeStep);
273 BondVector.Normalize();
274 if (BondVector.ScalarProduct(currentGradient) < 0) {
275 removed_bonds.push_back(*bonditer);
276#ifndef NDEBUG
277 const bool status =
278#endif
279 BGcreator.removeEdge(current_bond.leftatom->getId(), current_bond.rightatom->getId());
280 ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
281 }
282 }
283 BoostGraphHelpers::Nodeset_t bondside_set = NodeGatherer((*iter)->getId(), max_distance);
284 const BreadthFirstSearchGatherer::distance_map_t &distance_map = NodeGatherer.getDistances();
285 std::sort(bondside_set.begin(), bondside_set.end());
286 // re-add those edges
287 for (std::vector<bond::ptr>::const_iterator bonditer = removed_bonds.begin();
288 bonditer != removed_bonds.end(); ++bonditer)
289 BGcreator.addEdge((*bonditer)->leftatom->getId(), (*bonditer)->rightatom->getId());
290
291 // apply PositionUpdate to all nodes in the bondside_set
292 for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set.begin();
293 setiter != bondside_set.end(); ++setiter) {
294 const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
295 = distance_map.find(*setiter);
296 ASSERT( diter != distance_map.end(),
297 "ForceAnnealing() - could not find distance to an atom.");
298 const double factor = pow(damping_factor, diter->second);
299 LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
300 << factor << "*" << PositionUpdate);
301 if (GatheredUpdates.count((*setiter))) {
302 GatheredUpdates[(*setiter)] += factor*PositionUpdate;
303 } else {
304 GatheredUpdates.insert(
305 std::make_pair(
306 (*setiter),
307 factor*PositionUpdate) );
308 }
309 }
310
311 // extract largest components for showing progress of annealing
312 for(size_t i=0;i<NDIM;++i)
313 maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
314 }
315 // apply the gathered updates
316 for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
317 iter != GatheredUpdates.end(); ++iter) {
318 const atomId_t &atomid = iter->first;
319 const Vector &update = iter->second;
320 atom* const walker = World::getInstance().getAtom(AtomById(atomid));
321 ASSERT( walker != NULL,
322 "ForceAnnealing() - walker with id "+toString(atomid)+" has suddenly disappeared.");
323 LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
324 << ", namely " << *walker);
325 walker->setPosition( walker->getPosition() + update );
326 }
327 }
328
329 /** Reset function to unset static entities and artificial velocities.
330 *
331 */
332 void reset()
333 {
334 currentDeltat = 0.;
335 currentStep = 0;
336 }
337
338private:
339 //!> contains the current step in relation to maxsteps
340 static size_t currentStep;
341 //!> contains the maximum number of steps, determines initial and final step with currentStep
342 size_t maxSteps;
343 static double currentDeltat;
344 //!> minimum deltat for internal while loop (adaptive step width)
345 static double MinimumDeltat;
346 //!> contains the maximum bond graph distance up to which shifts of a single atom are spread
347 const int max_distance;
348 //!> the shifted is dampened by this factor with the power of the bond graph distance to the shift causing atom
349 const double damping_factor;
350};
351
352template <class T>
353double ForceAnnealing<T>::currentDeltat = 0.;
354template <class T>
355size_t ForceAnnealing<T>::currentStep = 0;
356template <class T>
357double ForceAnnealing<T>::MinimumDeltat = 1e-8;
358
359#endif /* FORCEANNEALING_HPP_ */
Note: See TracBrowser for help on using the repository browser.