240 {
241 long effectiveMaxIter = max_iterations_override;
242 QUILL_LOG_DEBUG(
243 log,
"Saddle point search started from reactant with energy {} eV.",
245
246 int optStatus;
247 bool firstIteration = true;
248 const char *forceLabel =
249 params.optimizer_options().convergence_metric_label.c_str();
250
251 if (
params.saddle_search_options().minmode_method ==
253 QUILL_LOG_DEBUG(
254 log,
"================= Using the GP Dimer Library =================");
257 QUILL_LOG_DEBUG(
log,
"GPR eigenvalue: {}",
261 }
263 QUILL_LOG_DEBUG(
log,
"[MinModeSaddleSearch] eigenvalue not negative");
265 }
267 params.saddle_search_options().zero_mode_abort_curvature) {
268 QUILL_LOG_DEBUG(
log,
"Zero mode eigenvalue: {}",
271 }
274 } else {
275
276 if (
params.saddle_search_options().minmode_method ==
279 "[Dimer] {:9s} {:9s} {:10s} {:18s} {:9s} "
280 "{:7s} {:6s} {:4s} {:5s}\n",
281 "Step", "Step Size", "Delta E", forceLabel, "Curvature",
282 "Torque", "Angle", "Rots", "Align");
283 }
else if (
params.saddle_search_options().minmode_method ==
285 QUILL_LOG_INFO(
287 "[Lanczos] {:9s} {:9s} {:10s} {:18s} {:9s} {:10s} {:7s} {:5s}\n",
288 "Step", "Step Size", "Delta E", forceLabel, "Curvature", "Rel Change",
289 "Angle", "Iters");
290 }
else if (
params.saddle_search_options().minmode_method ==
293 "[GPRDimer] {:9s} {:9s} {:10s} {:18s} {:9s} "
294 " {:7s} {:6s} {:4s}\n",
295 "Step", "Step Size", "Delta E", forceLabel, "Curvature",
296 "Torque", "Angle", "Rots");
297 }
298
299 std::string climbLabel = "climb";
300 std::string climbDatFilename = "climb.dat";
301
303
304 auto objf = std::make_shared<MinModeObjectiveFunction>(
306 auto write_climb_frame = [&](uint64_t frameIndex, bool append,
307 double stepSize, double de, double conv,
308 double eigenval, double torque, double angle,
309 long rotations) {
310 eonc::io::ConFrameMetadata metadata;
313 metadata.
scalars.push_back({
"step_size", stepSize});
314 metadata.
scalars.push_back({
"delta_e", de});
315 metadata.
scalars.push_back({
"convergence", conv});
316 metadata.
scalars.push_back({
"eigenvalue", eigenval});
317 metadata.
scalars.push_back({
"torque", torque});
318 metadata.
scalars.push_back({
"angle", angle});
319 metadata.
scalars.push_back({
"rotations",
static_cast<double>(rotations)});
322 }
323 if (
params.debug_options().write_movies) {
325 matter->matter2con(climbLabel, append, &metadata))) {
326 QUILL_LOG_WARNING(
log,
"Failed to write climb movie frame {}",
327 climbLabel);
328 }
332 }
333
334 if (
params.debug_options().write_deprecated_outs) {
335 std::ofstream climbDat(climbDatFilename,
336 append ? (std::ios::binary | std::ios::app)
337 : std::ios::binary);
338 if (climbDat) {
339 if (!append) {
340 climbDat << "iteration\tstep_size\tdelta_e\tconvergence"
341 "\teigenvalue\ttorque\tangle\trotations\n";
342 }
343 climbDat << std::format("{}\t{:.7e}\t{:.6f}\t{:.5e}\t{:.6f}"
344 "\t{:.6f}\t{:.4f}\t{}\n",
345 frameIndex, stepSize, de, conv, eigenval,
346 torque, angle, rotations);
347 }
348 }
349 };
351 write_climb_frame(0, false, 0.0, 0.0, objf->getConvergence(),
353 0);
354 }
355 if (
params.saddle_search_options().nonnegative_displacement_abort) {
356 objf->getGradient();
358 QUILL_LOG_DEBUG(
log,
"Nonnegative eigenvalue: {}",
362 }
363 }
364
367
368 while (!objf->isConverged() ||
iteration == 0) {
369
370 if (!firstIteration) {
371
372 if (
params.saddle_search_options().nonlocal_count_abort != 0) {
374 initialPosition -
matter->getPositions(),
375 params.saddle_search_options().nonlocal_distance_abort);
376 if (nm >=
params.saddle_search_options().nonlocal_count_abort) {
378 break;
379 }
380 }
381
383 params.saddle_search_options().zero_mode_abort_curvature) {
384 QUILL_LOG_DEBUG(
log,
"Zero mode eigenvalue: {}",
387 break;
388 }
389 }
390 firstIteration = false;
391
394 break;
395 }
396
398
399 try {
400 if (
params.saddle_search_options().confine_positive.bowl_breakout &&
403 optStatus = optim->step(-
params.optimizer_options().max_move);
404 } else {
405 optStatus = optim->step(
params.optimizer_options().max_move);
406 }
407 } catch (const eonc::DimerModeRestoredException &) {
408 QUILL_LOG_DEBUG(
409 log,
"Dimer restored to best state. Checking convergence...");
411 break;
412 } catch (const eonc::DimerModeLostException &) {
413 QUILL_LOG_WARNING(
log,
"Dimer lost mode completely. Aborting.");
415 break;
416 }
417
418 if (optStatus < 0) {
420 break;
421 }
422
424
425
426 if (
params.saddle_search_options().remove_rotation) {
428 }
429 double stepSize = (
matter->pbc(
matter->getPositions() - pos)).norm();
430
432
433
438 double conv = objf->getConvergence();
439
440 if (
params.saddle_search_options().minmode_method ==
442 QUILL_LOG_DEBUG(
444 "[Lanczos] {:9} {:9.6f} {:10.4f} {:18.5e} {:9.4f} {:10.6f} "
445 "{:7.3f} {:5}\n",
446 iteration, stepSize, de, conv, eigenval, torque, angle, rotations);
447 } else {
448 QUILL_LOG_DEBUG(
450 "[Dimer] {:9} {:9.7f} {:10.4f} {:18.5e} {:9.4f} {:7.3f} "
451 " {:6.3f} {:4}\n",
452 iteration, stepSize, de, conv, eigenval, torque, angle, rotations);
453 }
454
456 write_climb_frame(
static_cast<uint64_t
>(
iteration),
true, stepSize, de,
457 conv, eigenval, torque, angle, rotations);
458 }
459
460 if (
params.main_options().checkpoint) {
462 matter->matter2con(
"displacement_cp.con",
false))) {
463 QUILL_LOG_WARNING(
log,
"Failed to write displacement_cp.con");
464 }
467 }
468
469 if (de >
params.saddle_search_options().max_energy) {
471 break;
472 }
473
474
476 if (!dimer->rotationDidConverge) {
480 QUILL_LOG_DEBUG(
log,
"Dimer restored to valid state. C_tau={:.4f}",
481 dimer->getEigenvalue());
482 }
483 break;
484 }
485 }
486 }
487
490 }
491
492
493
494 const bool climbConverged = objf->isConverged();
495 const int statusBeforeGuard =
status;
498 QUILL_LOG_WARNING(
log,
"[MinModeSaddleSearch] objective not converged; "
499 "refusing STATUS_GOOD");
500 }
501
503 QUILL_LOG_DEBUG(
log,
"[MinModeSaddleSearch] eigenvalue not negative");
505 }
507 }
508
510}
Eigen::Matrix< double, Eigen::Dynamic, 3, eOnStorageOrder > AtomMatrix
static const char MINMODE_GPRDIMER[]
static const char MINMODE_DIMER[]
static const char MINMODE_LANCZOS[]
bool retain_climb_frames_
static constexpr int finalizeClimbStatus(int climbStatus, bool objectiveConverged) noexcept
Issue #20 policy: never leave climb as STATUS_GOOD when the climb objective is still unconverged (unf...
long numAtomsMoved(const AtomMatrix v1, double cutoff)
void rotationRemove(const AtomMatrix r1, std::shared_ptr< Matter > m2)
std::unique_ptr< Optimizer > mkOptim(std::shared_ptr< ObjectiveFunction > a_objf, OptType a_otype, const Parameters &a_params)
void saveMode(FILE *modeFile, std::shared_ptr< Matter > matter, AtomMatrix mode)
Write a mode; constrained axes are emitted as 0.
readcon::ConFrame matterToConFrame(Matter &m, const ConFrameMetadata *metadata)
Build a single stamped ConFrame from Matter (same builder as matter2con).
constexpr bool io_ok(IoStatus s) noexcept
void eigenmodeCompute(LowestEigenmode &s, std::shared_ptr< Matter > matter, AtomMatrix direction)
long eigenmodeTotalForceCalls(LowestEigenmode &s)
long eigenmodeStatsRotations(LowestEigenmode &s)
long eigenmodeTotalIterations(LowestEigenmode &s)
double eigenmodeStatsTorque(LowestEigenmode &s)
double eigenmodeStatsAngle(LowestEigenmode &s)