Changeset 83956e for src/Dynamics
- Timestamp:
- Jun 20, 2018, 8:20:42 AM (7 years ago)
- Branches:
- AutomationFragmentation_failures, Candidate_v1.6.1, ChemicalSpaceEvaluator, Exclude_Hydrogens_annealWithBondGraph, ForceAnnealing_with_BondGraph_contraction-expansion, StoppableMakroAction
- Children:
- 646f73
- Parents:
- 6a5921
- git-author:
- Frederik Heber <frederik.heber@…> (08/10/17 15:35:45)
- git-committer:
- Frederik Heber <frederik.heber@…> (06/20/18 08:20:42)
- Location:
- src/Dynamics
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
src/Dynamics/BondVectors.cpp
r6a5921 r83956e 285 285 } 286 286 287 void BondVectors::getProjectedGradientsForAtomAtStep( 288 const atom &_walker, 289 const Vector &_walkerGradient, 290 const size_t _timestep, 291 std::vector< std::vector<double> > &_projected_forces) const 292 { 293 // gather subset of BondVectors for the current atom 294 const std::vector<Vector> BondVectors = getAtomsBondVectorsAtStep(_walker, _timestep); 295 const BondList& ListOfBonds = _walker.getListOfBonds(); // we always use current bonds 296 297 // go over all its bonds 298 std::vector<Vector>::const_iterator vectoriter = BondVectors.begin(); 299 for(BondList::const_iterator bonditer = ListOfBonds.begin(); 300 bonditer != ListOfBonds.end(); ++bonditer, ++vectoriter) { 301 const bond::ptr ¤t_bond = *bonditer; 302 const Vector &BondVector = *vectoriter; 303 304 // bv goes from rightatom to leftatom 305 // as plus sign in force indicates expansion, minus indicates contraction 306 // leftatom: from right to left means same sign expansion, opposite contraction 307 // rightatom: from right to left means opposite sign expansion, same sign contraction 308 const double sign = (&_walker == current_bond->leftatom) ? 1. : -1.; 309 const double temp = sign*_walkerGradient.ScalarProduct(BondVector); 310 LOG(4, "DEBUG: BondVector " << BondVector << " receives projected force of " 311 << sign << "*" << _walkerGradient << "*" << BondVector << " = " << temp); 312 const size_t index = getIndexForBond(current_bond); 313 ASSERT( index != (size_t)-1, 314 "ForceAnnealing() - could not find bond "+toString(*current_bond) 315 +" in bondvectors"); 316 std::vector<double> &forcelist = (&_walker == current_bond->leftatom) ? 317 _projected_forces[leftside] : _projected_forces[rightside]; 318 forcelist[index] = temp; 319 } 320 ASSERT( vectoriter == BondVectors.end(), 321 "BondVectors::getRemnantGradientForAtomAtStep() - vectoriter is not at end when it should be."); 322 } 323 287 324 Vector BondVectors::getRemnantGradientForAtomAtStep( 288 325 const atom &_walker, -
src/Dynamics/BondVectors.hpp
r6a5921 r83956e 56 56 const size_t &_step); 57 57 58 /** Returns the number of bonds whose bond vectors have been registered. 59 * 60 * \param number of bonds 61 */ 62 size_t size() const 63 { return container.size(); } 64 58 65 /** Getter for the sorted bonds. 59 66 * … … 144 151 forcestore_t _forcestore) const; 145 152 153 /** Calculates the \a _walkkerGradient projected onto the bond vector for 154 * every of \a _walker's bonds. 155 * 156 * \param _walker atom 157 * \param _walkerGradient atom's gradient 158 * \param _timestep time step 159 * \param _projected_forces list of forces for every \a whichatom_t and every bond 160 */ 161 void getProjectedGradientsForAtomAtStep( 162 const atom &_walker, 163 const Vector &_walkerGradient, 164 const size_t _timestep, 165 std::vector< std::vector<double> > &_projected_forces) const; 166 167 // knowing the number of bonds in total, we can setup the storage for the 168 // projected forces 169 enum whichatom_t { 170 leftside=0, 171 rightside=1, 172 MAX_sides 173 }; 174 146 175 private: 147 176 /** Calculates the bond vector for each bond in the internal container. … … 177 206 //!> internal map for bond Bondvector association 178 207 mutable mapped_t current_mapped_vectors; 208 179 209 }; 180 210 -
src/Dynamics/ForceAnnealing.hpp
r6a5921 r83956e 207 207 "ForceAnnealing::anneal() - if currentStep is "+toString(currentStep) 208 208 +", then there should be at least three time steps."); 209 const Vector oldPosition = (*iter)->getPositionAtStep(OldTimeStep);209 const Vector &oldPosition = (*iter)->getPositionAtStep(OldTimeStep); 210 210 const Vector PositionDifference = currentPosition - oldPosition; 211 211 LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition); … … 259 259 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) { 260 260 // atom's force vector gives steepest descent direction 261 const Vector oldPosition = (*iter)->getPositionAtStep(OldTimeStep);262 const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);263 const Vector oldGradient = (*iter)->getAtomicForceAtStep(OldTimeStep);264 const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);261 const Vector &oldPosition = (*iter)->getPositionAtStep(OldTimeStep); 262 const Vector ¤tPosition = (*iter)->getPositionAtStep(CurrentTimeStep); 263 const Vector &oldGradient = (*iter)->getAtomicForceAtStep(OldTimeStep); 264 const Vector ¤tGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep); 265 265 LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition); 266 266 LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition); … … 304 304 } 305 305 306 // knowing the number of bonds in total, we can setup the storage for the307 // projected forces308 enum whichatom_t {309 leftside=0,310 rightside=1,311 MAX_sides312 };313 314 /** Helper function to put bond force into a container.315 *316 * \param _walker atom317 * \param _current_bond current bond of \a _walker318 * \param _timestep time step319 * \param _force calculated bond force320 * \param _bv bondvectors for obtaining the correct index321 * \param _projected_forces container322 */323 static void ForceStore(324 const atom &_walker,325 const bond::ptr &_current_bond,326 const size_t &_timestep,327 const double _force,328 const BondVectors &_bv,329 std::vector< // time step330 std::vector< // which bond side331 std::vector<double> > // over all bonds332 > &_projected_forces)333 {334 std::vector<double> &forcelist = (&_walker == _current_bond->leftatom) ?335 _projected_forces[_timestep][leftside] : _projected_forces[_timestep][rightside];336 const size_t index = _bv.getIndexForBond(_current_bond);337 ASSERT( index != (size_t)-1,338 "ForceAnnealing() - could not find bond "+toString(*_current_bond)339 +" in bondvectors");340 forcelist[index] = _force;341 }342 343 306 /** Performs Gradient optimization on the bonds with BarzilaiBorwein stepwdith. 344 307 * … … 377 340 BreadthFirstSearchGatherer NodeGatherer(BGcreator); 378 341 379 /// We assume that a force is local, i.e. a bond is too short yet and hence 380 /// the atom needs to be moved. However, all the adjacent (bound) atoms might 381 /// already be at the perfect distance. If we just move the atom alone, we ruin 382 /// all the other bonds. Hence, it would be sensible to move every atom found 383 /// through the bond graph in the direction of the force as well by the same 384 /// PositionUpdate. This is almost what we are going to do. 385 386 /// One issue is: If we need to shorten bond, then we use the PositionUpdate 387 /// also on the the other bond partner already. This is because it is in the 388 /// direction of the bond. Therefore, the update is actually performed twice on 389 /// each bond partner, i.e. the step size is twice as large as it should be. 390 /// This problem only occurs when bonds need to be shortened, not when they 391 /// need to be made longer (then the force vector is facing the other 392 /// direction than the bond vector). 393 /// As a remedy we need to average the force on either end of the bond and 394 /// check whether each gradient points inwards out or outwards with respect 395 /// to the bond and then shift accordingly. 396 397 /// One more issue is that the projection onto the bond directions does not 398 /// recover the gradient but may be larger as the bond directions are a 399 /// generating system and not a basis (e.g. 3 bonds on a plane where 2 would 400 /// suffice to span the plane). To this end, we need to account for the 401 /// overestimation and obtain a weighting for each bond. 342 /** We assume that a force is local, i.e. a bond is too short yet and hence 343 * the atom needs to be moved. However, all the adjacent (bound) atoms might 344 * already be at the perfect distance. If we just move the atom alone, we ruin 345 * all the other bonds. Hence, it would be sensible to move every atom found 346 * through the bond graph in the direction of the force as well by the same 347 * PositionUpdate. This is almost what we are going to do, see below. 348 * 349 * This is to make the force a little more global in the sense of a multigrid 350 * solver that uses various coarser grids to transport errors more effectively 351 * over finely resolved grids. 352 * 353 */ 354 355 /** The idea is that we project the gradients onto the bond vectors and determine 356 * from the sum of projected gradients from either side whether the bond is 357 * to contract or to expand. As the gradient acting as the normal vector of 358 * a plane supported at the position of the atom separates all bonds into two 359 * sets, we check whether all on one side are contracting and all on the other 360 * side are expanding. In this case we may move not only the atom itself but 361 * may propagate its update along a limited-horizon BFS to neighboring atoms. 362 * 363 */ 402 364 403 365 // initialize helper class for bond vectors using bonds from range of atoms … … 407 369 AtomicForceManipulator<T>::atoms.end(), 408 370 _TimeStep); // use time step to update here as this is the current set of bonds 409 const BondVectors::container_t &sorted_bonds = bv.getSorted(); 410 411 std::vector< // time step 412 std::vector< // which bond side 413 std::vector<double> > // over all bonds 414 > projected_forces(2); // one for leftatoms, one for rightatoms (and for both time steps) 415 for (size_t i=0;i<2;++i) { 416 projected_forces[i].resize(MAX_sides); 417 for (size_t j=0;j<MAX_sides;++j) 418 projected_forces[i][j].resize(sorted_bonds.size(), 0.); 419 } 420 421 // for each atom we need to gather weights and then project the gradient 422 typedef std::map<atomId_t, BondVectors::weights_t > weights_per_atom_t; 423 std::vector<weights_per_atom_t> weights_per_atom(2); 424 typedef std::map<atomId_t, Vector> RemnantGradient_per_atom_t; 425 RemnantGradient_per_atom_t RemnantGradient_per_atom; 426 for (size_t timestep = 0; timestep <= 1; ++timestep) { 427 const size_t ReferenceTimeStep = CurrentTimeStep-timestep; 428 LOG(2, "DEBUG: given time step is " << _TimeStep 429 << ", timestep is " << timestep 430 << ", and ReferenceTimeStep is " << ReferenceTimeStep); 431 432 for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin(); 433 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) { 434 const atom &walker = *(*iter); 435 const Vector &walkerGradient = walker.getAtomicForceAtStep(ReferenceTimeStep); 436 LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely " 437 << walker << " is " << walkerGradient << " with magnitude of " 438 << walkerGradient.Norm()); 439 440 const BondList& ListOfBonds = walker.getListOfBonds(); 441 if (walkerGradient.Norm() > MYEPSILON) { 442 443 // gather subset of BondVectors for the current atom 444 const std::vector<Vector> BondVectors = 445 bv.getAtomsBondVectorsAtStep(walker, ReferenceTimeStep); 446 447 // go through all its bonds and calculate what magnitude is represented 448 // by the others i.e. sum of scalar products against other bonds 449 const std::pair<weights_per_atom_t::iterator, bool> inserter = 450 weights_per_atom[timestep].insert( 451 std::make_pair(walker.getId(), 452 bv.getWeightsForAtomAtStep(walker, BondVectors, ReferenceTimeStep)) ); 453 ASSERT( inserter.second, 454 "ForceAnnealing::operator() - weight map for atom "+toString(walker) 455 +" and time step "+toString(timestep)+" already filled?"); 456 BondVectors::weights_t &weights = inserter.first->second; 457 ASSERT( weights.size() == ListOfBonds.size(), 458 "ForceAnnealing::operator() - number of weights " 459 +toString(weights.size())+" does not match number of bonds " 460 +toString(ListOfBonds.size())+", error in calculation?"); 461 462 // projected gradient over all bonds and place in one of projected_forces 463 // using the obtained weights 464 BondVectors::forcestore_t forcestoring = 465 boost::bind(&ForceAnnealing::ForceStore, _1, _2, _3, _4, 466 boost::cref(bv), boost::ref(projected_forces)); 467 const Vector RemnantGradient = bv.getRemnantGradientForAtomAtStep( 468 walker, walkerGradient, BondVectors, weights, timestep, forcestoring 469 ); 470 RemnantGradient_per_atom.insert( std::make_pair(walker.getId(), RemnantGradient) ); 471 } else { 472 LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than " 473 << MYEPSILON << " for atom " << walker); 474 // note that projected_forces is initialized to full length and filled 475 // with zeros. Hence, nothing to do here 371 372 std::vector< // which bond side 373 std::vector<double> > // over all bonds 374 projected_forces; // one for leftatoms, one for rightatoms 375 projected_forces.resize(BondVectors::MAX_sides); 376 for (size_t j=0;j<BondVectors::MAX_sides;++j) 377 projected_forces[j].resize(bv.size(), 0.); 378 379 // for each atom we need to project the gradient 380 for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin(); 381 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) { 382 const atom &walker = *(*iter); 383 const Vector &walkerGradient = walker.getAtomicForceAtStep(CurrentTimeStep); 384 const double GradientNorm = walkerGradient.Norm(); 385 LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely " 386 << walker << " is " << walkerGradient << " with magnitude of " 387 << GradientNorm); 388 389 if (GradientNorm > MYEPSILON) { 390 bv.getProjectedGradientsForAtomAtStep( 391 walker, walkerGradient, CurrentTimeStep, projected_forces 392 ); 393 } else { 394 LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than " 395 << MYEPSILON << " for atom " << walker); 396 // note that projected_forces is initialized to full length and filled 397 // with zeros. Hence, nothing to do here 398 } 399 } 400 401 std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end 402 for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin(); 403 iter != AtomicForceManipulator<T>::atoms.end(); ++iter) { 404 atom &walker = *(*iter); 405 406 /// calculate step width 407 const Vector &oldPosition = (*iter)->getPositionAtStep(OldTimeStep); 408 const Vector ¤tPosition = (*iter)->getPositionAtStep(CurrentTimeStep); 409 const Vector &oldGradient = (*iter)->getAtomicForceAtStep(OldTimeStep); 410 const Vector ¤tGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep); 411 LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition); 412 LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition); 413 LOG(4, "DEBUG: oldGradient for atom #" << (*iter)->getId() << " is " << oldGradient); 414 LOG(4, "DEBUG: currentGradient for atom #" << (*iter)->getId() << " is " << currentGradient); 415 // LOG(4, "DEBUG: Force for atom #" << (*iter)->getId() << " is " << currentGradient); 416 417 // we use Barzilai-Borwein update with position reversed to get descent 418 const Vector PositionDifference = currentPosition - oldPosition; 419 const Vector GradientDifference = (currentGradient - oldGradient); 420 double stepwidth = 0.; 421 if (GradientDifference.Norm() > MYEPSILON) 422 stepwidth = fabs(PositionDifference.ScalarProduct(GradientDifference))/ 423 GradientDifference.NormSquared(); 424 if (fabs(stepwidth) < 1e-10) { 425 // dont' warn in first step, deltat usage normal 426 if (currentStep != 1) 427 ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead."); 428 stepwidth = currentDeltat; 429 } 430 Vector PositionUpdate = stepwidth * currentGradient; 431 LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate); 432 433 /** go through all bonds and check projected_forces and side of plane 434 * the idea is that if all bonds on one side are contracting ones or expanding, 435 * respectively, then we may shift not only the atom with respect to its 436 * gradient but also its neighbors (towards contraction or towards 437 * expansion depending on direction of gradient). 438 * if they are mixed on both sides of the plane, then we simply shift 439 * only the atom itself. 440 * if they are not mixed on either side, then we also only shift the 441 * atom, namely away from expanding and towards contracting bonds. 442 */ 443 444 // sign encodes side of plane and also encodes contracting(-) or expanding(+) 445 typedef std::vector<int> sides_t; 446 typedef std::vector<int> types_t; 447 sides_t sides; 448 types_t types; 449 const BondList& ListOfBonds = walker.getListOfBonds(); 450 for(BondList::const_iterator bonditer = ListOfBonds.begin(); 451 bonditer != ListOfBonds.end(); ++bonditer) { 452 const bond::ptr ¤t_bond = *bonditer; 453 454 // BondVector goes from bond::rightatom to bond::leftatom 455 const size_t index = bv.getIndexForBond(current_bond); 456 std::vector<double> &forcelist = (&walker == current_bond->leftatom) ? 457 projected_forces[BondVectors::leftside] : projected_forces[BondVectors::rightside]; 458 const double &temp = forcelist[index]; 459 sides.push_back( -1.*temp/fabs(temp) ); // BondVectors has exactly opposite sign for sides decision 460 const double sum = 461 projected_forces[BondVectors::leftside][index]+projected_forces[BondVectors::rightside][index]; 462 types.push_back( sum/fabs(sum) ); 463 LOG(4, "DEBUG: Bond " << *current_bond << " is on side " << sides.back() 464 << " and has type " << types.back()); 465 } 466 // /// check whether both conditions are compatible: 467 // // i.e. either we have ++/-- for all entries in sides and types 468 // // or we have +-/-+ for all entries 469 // // hence, multiplying and taking the sum and its absolute value 470 // // should be equal to the maximum number of entries 471 // sides_t results; 472 // std::transform( 473 // sides.begin(), sides.end(), 474 // types.begin(), 475 // std::back_inserter(results), 476 // std::multiplies<int>); 477 // int result = abs(std::accumulate(results.begin(), results.end(), 0, std::plus<int>)); 478 479 std::vector<size_t> first_per_side(2, (size_t)-1); //!< mark down one representative from either side 480 std::vector< std::vector<int> > types_per_side(2); //!< gather all types on each side 481 types_t::const_iterator typesiter = types.begin(); 482 for (sides_t::const_iterator sidesiter = sides.begin(); 483 sidesiter != sides.end(); ++sidesiter, ++typesiter) { 484 const size_t index = (*sidesiter+1)/2; 485 types_per_side[index].push_back(*typesiter); 486 if (first_per_side[index] == (size_t)-1) 487 first_per_side[index] = std::distance(const_cast<const sides_t &>(sides).begin(), sidesiter); 488 } 489 LOG(4, "DEBUG: First on side minus is " << first_per_side[0] << ", and first on side plus is " 490 << first_per_side[1]); 491 //!> enumerate types per side with a little witching with the numbers to allow easy setting from types 492 enum whichtypes_t { 493 contracting=0, 494 unset=1, 495 expanding=2, 496 mixed 497 }; 498 std::vector<int> typeside(2, unset); 499 for(size_t i=0;i<2;++i) { 500 for (std::vector<int>::const_iterator tpsiter = types_per_side[i].begin(); 501 tpsiter != types_per_side[i].end(); ++tpsiter) { 502 if (typeside[i] == unset) { 503 typeside[i] = *tpsiter+1; //contracting(0) or expanding(2) 504 } else { 505 if (typeside[i] != (*tpsiter+1)) // no longer he same type 506 typeside[i] = mixed; 507 } 476 508 } 477 509 } 478 } 479 480 // step through each bond and shift the atoms 481 std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end 482 483 LOG(3, "DEBUG: current step is " << currentStep << ", given time step is " << CurrentTimeStep); 484 const BondVectors::mapped_t bondvectors = bv.getBondVectorsAtStep(CurrentTimeStep); 485 486 for (BondVectors::container_t::const_iterator bondsiter = sorted_bonds.begin(); 487 bondsiter != sorted_bonds.end(); ++bondsiter) { 488 const bond::ptr ¤t_bond = *bondsiter; 489 const size_t index = bv.getIndexForBond(current_bond); 490 const atom* bondatom[MAX_sides] = { 491 current_bond->leftatom, 492 current_bond->rightatom 493 }; 494 495 // remove the edge 510 LOG(4, "DEBUG: Minus side is " << typeside[0] << " and plus side is " << typeside[1]); 511 512 typedef std::vector< std::pair<atomId_t, atomId_t> > RemovedEdges_t; 513 if ((typeside[0] != mixed) || (typeside[1] != mixed)) { 514 const size_t sideno = ((typeside[0] != mixed) && (typeside[0] != unset)) ? 0 : 1; 515 LOG(4, "DEBUG: Chosen side is " << sideno << " with type " << typeside[sideno]); 516 ASSERT( (typeside[sideno] == contracting) || (typeside[sideno] == expanding), 517 "annealWithBondGraph_BB() - chosen side is either expanding nor contracting."); 518 // one side is not mixed, all bonds on one side are of same type 519 // hence, find out which bonds to exclude 520 const BondList& ListOfBonds = walker.getListOfBonds(); 521 522 // away from gradient (minus) and contracting 523 // or towards gradient (plus) and expanding 524 // gather all on same side and remove 525 const double sign = 526 (sides[first_per_side[sideno]] == types[first_per_side[sideno]]) 527 ? sides[first_per_side[sideno]] : -1.*sides[first_per_side[sideno]]; 528 LOG(4, "DEBUG: Removing edges from side with sign " << sign); 529 BondList::const_iterator bonditer = ListOfBonds.begin(); 530 RemovedEdges_t RemovedEdges; 531 for (sides_t::const_iterator sidesiter = sides.begin(); 532 sidesiter != sides.end(); ++sidesiter, ++bonditer) { 533 if (*sidesiter == sign) { 534 // remove the edge 535 const bond::ptr ¤t_bond = *bonditer; 536 LOG(5, "DEBUG: Removing edge " << *current_bond); 537 RemovedEdges.push_back( std::make_pair( 538 current_bond->leftatom->getId(), 539 current_bond->rightatom->getId()) 540 ); 496 541 #ifndef NDEBUG 497 const bool status =542 const bool status = 498 543 #endif 499 BGcreator.removeEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId()); 500 ASSERT( status, "ForceAnnealing() - edge to found bond is not present?"); 501 502 // gather nodes for either atom 503 BoostGraphHelpers::Nodeset_t bondside_set[MAX_sides]; 504 BreadthFirstSearchGatherer::distance_map_t distance_map[MAX_sides]; 505 for (size_t side=leftside;side<MAX_sides;++side) { 506 bondside_set[side] = NodeGatherer(bondatom[side]->getId(), max_distance); 507 distance_map[side] = NodeGatherer.getDistances(); 508 std::sort(bondside_set[side].begin(), bondside_set[side].end()); 509 } 510 511 // re-add edge 512 BGcreator.addEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId()); 513 514 // do for both leftatom and rightatom of bond 515 for (size_t side = leftside; side < MAX_sides; ++side) { 516 const double &bondforce = projected_forces[0][side][index]; 517 const double &oldbondforce = projected_forces[1][side][index]; 518 const double bondforcedifference = fabs(bondforce - oldbondforce); 519 LOG(4, "DEBUG: bondforce for " << (side == leftside ? "left" : "right") 520 << " side of bond is " << bondforce); 521 LOG(4, "DEBUG: oldbondforce for " << (side == leftside ? "left" : "right") 522 << " side of bond is " << oldbondforce); 523 // if difference or bondforce itself is zero, do nothing 524 if ((fabs(bondforce) < MYEPSILON) || (fabs(bondforcedifference) < MYEPSILON)) 525 continue; 526 527 // get BondVector to bond 528 const BondVectors::mapped_t::const_iterator bviter = 529 bondvectors.find(current_bond); 530 ASSERT( bviter != bondvectors.end(), 531 "ForceAnnealing() - cannot find current_bond ?"); 532 ASSERT( fabs(bviter->second.Norm() -1.) < MYEPSILON, 533 "ForceAnnealing() - norm of BondVector is not one"); 534 const Vector &BondVector = bviter->second; 535 536 // calculate gradient and position differences for stepwidth 537 const Vector currentGradient = bondforce * BondVector; 538 LOG(4, "DEBUG: current projected gradient for " 539 << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient); 540 const Vector &oldPosition = bondatom[side]->getPositionAtStep(OldTimeStep); 541 const Vector ¤tPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep); 542 const Vector PositionDifference = currentPosition - oldPosition; 543 LOG(4, "DEBUG: old position is " << oldPosition); 544 LOG(4, "DEBUG: current position is " << currentPosition); 545 LOG(4, "DEBUG: difference in position is " << PositionDifference); 546 LOG(4, "DEBUG: bondvector is " << BondVector); 547 const double projected_PositionDifference = PositionDifference.ScalarProduct(BondVector); 548 LOG(4, "DEBUG: difference in position projected onto bondvector is " 549 << projected_PositionDifference); 550 LOG(4, "DEBUG: abs. difference in forces is " << bondforcedifference); 551 552 // calculate step width 553 double stepwidth = 554 fabs(projected_PositionDifference)/bondforcedifference; 555 if (fabs(stepwidth) < 1e-10) { 556 // dont' warn in first step, deltat usage normal 557 if (currentStep != 1) 558 ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead."); 559 stepwidth = currentDeltat; 544 BGcreator.removeEdge(RemovedEdges.back()); 545 ASSERT( status, "ForceAnnealing() - edge to found bond is not present?"); 546 } 560 547 } 561 Vector PositionUpdate = stepwidth * currentGradient; 562 LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate); 563 564 // add PositionUpdate for all nodes in the bondside_set 565 for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set[side].begin(); 566 setiter != bondside_set[side].end(); ++setiter) { 548 // perform limited-horizon BFS 549 BoostGraphHelpers::Nodeset_t bondside_set; 550 BreadthFirstSearchGatherer::distance_map_t distance_map; 551 bondside_set = NodeGatherer(walker.getId(), max_distance); 552 distance_map = NodeGatherer.getDistances(); 553 std::sort(bondside_set.begin(), bondside_set.end()); 554 555 // re-add edge 556 for (RemovedEdges_t::const_iterator edgeiter = RemovedEdges.begin(); 557 edgeiter != RemovedEdges.end(); ++edgeiter) 558 BGcreator.addEdge(edgeiter->first, edgeiter->second); 559 560 // update position with dampening factor on the discovered bonds 561 for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set.begin(); 562 setiter != bondside_set.end(); ++setiter) { 567 563 const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter 568 = distance_map [side].find(*setiter);569 ASSERT( diter != distance_map [side].end(),564 = distance_map.find(*setiter); 565 ASSERT( diter != distance_map.end(), 570 566 "ForceAnnealing() - could not find distance to an atom."); 571 567 const double factor = pow(damping_factor, diter->second+1); … … 581 577 } 582 578 } 579 } else { 580 // simple atomic annealing, i.e. damping factor of 1 581 LOG(3, "DEBUG: Update for atom #" << walker.getId() << " will be " << PositionUpdate); 582 GatheredUpdates.insert( 583 std::make_pair( 584 walker.getId(), 585 PositionUpdate) ); 583 586 } 584 587 } … … 593 596 } 594 597 595 // remove center of weight translation from gathered updates596 Vector CommonTranslation;597 for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();598 iter != GatheredUpdates.end(); ++iter) {599 const Vector &update = iter->second;600 CommonTranslation += update;601 }602 CommonTranslation *= 1./(double)GatheredUpdates.size();603 LOG(3, "DEBUG: Subtracting common translation " << CommonTranslation604 << " from all updates.");598 // // remove center of weight translation from gathered updates 599 // Vector CommonTranslation; 600 // for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin(); 601 // iter != GatheredUpdates.end(); ++iter) { 602 // const Vector &update = iter->second; 603 // CommonTranslation += update; 604 // } 605 // CommonTranslation *= 1./(double)GatheredUpdates.size(); 606 // LOG(3, "DEBUG: Subtracting common translation " << CommonTranslation 607 // << " from all updates."); 605 608 606 609 // apply the gathered updates and set remnant gradients for atomic annealing … … 615 618 << ", namely " << *walker); 616 619 walker->setPositionAtStep(_TimeStep, 617 walker->getPositionAtStep(CurrentTimeStep) 618 + update - CommonTranslation); 619 // walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] ); 620 walker->getPositionAtStep(CurrentTimeStep) + update); // - CommonTranslation); 620 621 } 621 622
Note:
See TracChangeset
for help on using the changeset viewer.