70 {
71
72 VectorXd initialDirection = VectorXd::Map(initialDirectionAtomMatrix.data(),
73 3 * matter->numberOfAtoms());
74 tau = initialDirection.array() * matter->getFreeV().array();
77 if (
tau.norm() > 1e-10) {
78 eonc::safemath::safe_normalize_inplace(
tau);
79 } else {
80
82 tau =
tau.array() * matter->getFreeV().array();
83 eonc::safemath::safe_normalize_inplace(
tau);
84 }
85
86
87 double bestNegativeCurvature = std::numeric_limits<double>::max();
88 VectorXd bestTau =
tau;
89 VectorXd bestX0Positions;
90 VectorXd bestG0, bestG1;
91
92
94
95 {
96
97 auto x1Pot =
x1->getPotential();
100 x1->setPotential(x1Pot);
101 }
102 VectorXd x0_r =
x0->getPositionsV();
103 bestX0Positions = x0_r;
104
105 double delta =
params.main_options().finiteDifference;
106 x1->setPositionsV(x0_r + delta *
tau);
107
108
109
110 {
111 Matter *
const ends[] = {
x0.get(),
x1.get()};
113 }
114
115
116 if (
x1->getPotentialEnergy() -
x0->getPotentialEnergy() > 10.0 * delta) {
118 x1->setPositionsV(x0_r + delta *
tau);
119 QUILL_LOG_DEBUG(
120 log,
"[IDimer] Initial tangent flipped due to high energy wall.");
121 }
122
123
126 AtomMatrix::Map(
tau.data(), matter->numberOfAtoms(), 3),
127 static_cast<quill::Logger *
>(
log))) {
128 C_tau = alt->eigenvalue;
129 tau = VectorXd::Map(alt->eigenvector.data(), 3 * matter->numberOfAtoms());
132 tau =
tau.array() * matter->getFreeV().array();
133 eonc::safemath::safe_normalize_inplace(
tau);
134 x0_r = matter->getPositionsV();
135 x0->setPositionsV(x0_r);
136 x1->setPositionsV(x0_r + delta *
tau);
140 return;
141 }
142
148 }
151 }
152
153 VectorXd x1_rp, x1_r, tau_prime, tau_Old, g1_prime;
154 double phi_tol =
156 double phi_prime = 0.0;
157 double phi_min = 0.0;
158
160
161
162 auto x1p = std::make_shared<Matter>(
x1->getPotential(),
params);
163
164
165 if (
params.dimer_options().remove_rotation) {
167 AtomMatrix::Map(x0_r.data(),
x0->numberOfAtoms(), 3),
x1);
168 x1_r =
x1->getPositionsV();
169 tau =
x1->pbcV(x1_r - x0_r);
170 eonc::safemath::safe_normalize_inplace(
tau);
171 x1_r = x0_r +
tau * delta;
172 x1->setPositionsV(x1_r);
173 }
174
175
176
177
178
179
180 VectorXd g0, g1;
181
182
184 (
pot->needsPerImageInstance() &&
185 x0->getPotential().get() !=
x1->getPotential().get());
186 if (
pot->supportsBatchEvaluation()) {
187 long n =
x0->numberOfAtoms();
188 bool x0dirty =
x0->needsForceUpdate();
189 bool x1dirty =
x1->needsForceUpdate();
190
191 if (x0dirty && x1dirty) {
192 auto nrs0 =
x0->getAtomicNrs();
193 auto nrs1 =
x1->getAtomicNrs();
194 auto box0 =
x0->getCell();
195 auto box1 =
x1->getCell();
196 const double *posVec[] = {
x0->getPositions().data(),
197 x1->getPositions().data()};
198 const int *nrsVec[] = {nrs0.data(), nrs1.data()};
199 double *frcVec[] = {
x0->forcesData(),
x1->forcesData()};
200 double energies[2], vars[2];
201 const double *boxVec[] = {box0.data(), box1.data()};
202 pot->forceBatch(2, n, posVec, nrsVec, frcVec, energies, vars, boxVec);
203 x0->setComputedPotential(energies[0], vars[0]);
204 x1->setComputedPotential(energies[1], vars[1]);
205 } else if (x1dirty) {
206 auto nrs =
x1->getAtomicNrs();
207 auto box =
x1->getCell();
208 const double *posVec[] = {
x1->getPositions().data()};
209 const int *nrsVec[] = {nrs.data()};
210 double *frcVec[] = {
x1->forcesData()};
211 double energies[1], vars[1];
212 const double *boxVec[] = {box.data()};
213 pot->forceBatch(1, n, posVec, nrsVec, frcVec, energies, vars, boxVec);
214 x1->setComputedPotential(energies[0], vars[0]);
215 } else if (x0dirty) {
217 }
218 g0 = -
x0->getForcesV();
219 g1 = -
x1->getForcesV();
220 }
else if (
params.main_options().parallel && canParallel) {
221
222
223
224
225 std::exception_ptr t0Error;
226 std::thread t0([&] {
227 try {
228 g0 = -
x0->getForcesV();
229 } catch (...) {
230 t0Error = std::current_exception();
231 }
232 });
233 try {
234 g1 = -
x1->getForcesV();
235 } catch (...) {
236 t0.join();
237 throw;
238 }
239 t0.join();
240 if (t0Error)
241 std::rethrow_exception(t0Error);
242 } else {
243 g0 = -
x0->getForcesV();
244 g1 = -
x1->getForcesV();
245 }
246
247 bestG0 = g0;
248 bestG1 = g1;
249
256
257 do {
258
259
260 F_R = -2.0 * (g1 - g0) + 2.0 * ((g1 - g0).dot(
tau)) *
tau;
262
263
265 theta = eonc::safemath::safe_normalized(
F_R);
266
271 } else {
276 : 0.0;
277 }
278
284 }
285 eonc::safemath::safe_normalize_inplace(
theta);
287
290 VectorXd s0 =
tau - tau_Old;
292 VectorXd y0 =
296 } else {
298 }
299
300 double H0 = 1.0 / 60.0;
301 size_t loopmax =
s.size();
302 std::vector<double> alpha(loopmax);
303
305 for (long i = static_cast<long>(loopmax) - 1; i >= 0; i--) {
306 alpha[i] =
rho[i] *
s[i].dot(q);
307 q -= alpha[i] *
y[i];
308 }
309 VectorXd z = H0 * q;
310 for (size_t i = 0; i < loopmax; i++) {
311 double bv =
rho[i] *
y[i].dot(z);
312 z +=
s[i] * (alpha[i] - bv);
313 }
314
315 double vd = std::clamp(-eonc::safemath::safe_normalized(z).dot(
316 eonc::safemath::safe_normalized(
F_R)),
317 -1.0, 1.0);
318 double angle =
320
321 if (angle > 87.0) {
326 }
327
328 theta = -eonc::safemath::safe_normalized(z);
330 eonc::safemath::safe_normalize_inplace(
theta);
331
335 }
336
337
339
340
341 if (
C_tau < bestNegativeCurvature) {
342 bestNegativeCurvature =
C_tau;
344 bestX0Positions =
x0->getPositionsV();
345 bestG0 = g0;
346 bestG1 = g1;
348 }
349
350
351 double d_C_tau_d_phi =
354 d_C_tau_d_phi, 2.0 * std::abs(
C_tau), 0.0);
356
357 double alignment = std::abs(
tau.dot(referenceMode));
358
359 if (std::abs(phi_prime) > phi_tol) {
360 double b1 = 0.5 * d_C_tau_d_phi;
361
362
363 x0_r =
x0->getPositionsV();
364 tau_prime =
tau * std::cos(phi_prime) +
theta * std::sin(phi_prime);
365 tau_prime = eonc::safemath::safe_normalized(tau_prime);
366 x1_rp = x0_r + tau_prime * delta;
367
369 x1p->setPositionsV(x1_rp);
370 g1_prime = -x1p->getForcesV();
371
374
375 double C_tau_prime =
377
378
380 C_tau - C_tau_prime + b1 * std::sin(2.0 * phi_prime),
381 1.0 - std::cos(2.0 * phi_prime), 0.0);
382 double a0 = 2.0 * (
C_tau - a1);
384
385 double C_tau_min = 0.5 * a0 + a1 * std::cos(2.0 * phi_min) +
386 b1 * std::sin(2.0 * phi_min);
387
388
389 if (C_tau_min >
C_tau) {
391 C_tau_min = 0.5 * a0 + a1 * std::cos(2.0 * phi_min) +
392 b1 * std::sin(2.0 * phi_min);
393 }
394
395
398 }
400
401
402 tau =
tau * std::cos(phi_min) +
theta * std::sin(phi_min);
403 tau = eonc::safemath::safe_normalized(
tau);
404 x1_r = x0_r +
tau * delta;
405
406
407 if (
params.dimer_options().remove_rotation) {
408 x1->setPositionsV(x1_r);
410 AtomMatrix::Map(x0_r.data(),
x0->numberOfAtoms(), 3),
x1);
411 x1_r =
x1->getPositionsV();
413 eonc::safemath::safe_normalize_inplace(
tau);
414 x1_r = x0_r +
tau * delta;
415 }
416
417 x1->setPositionsV(x1_r);
419
420
421 double sin_pp = std::sin(phi_prime);
423 0.0) +
425 g0 * (1.0 - std::cos(phi_min) -
426 std::sin(phi_min) * std::tan(phi_prime * 0.5));
427
430 QUILL_LOG_INFO(
432 "[IDimerRot] ----- --------- ---------- ------------------ "
433 "{:9.4f} {:7.3f} {:6.3f} {:4} {:5.3f}",
435 } else {
436 QUILL_LOG_INFO(
438 "[IDimerRot] ----- --------- ---------- ------------------ "
439 "{:9.4f} {:7.3f} ------ ---- {:5.3f}",
440 C_tau,
F_R.norm() / delta, alignment);
441 }
442
443
444 if (alignment <
params.neb_options().climbing_image.ocineb.angle_tol &&
445 params.neb_options().climbing_image.ocineb.use_mmf) {
446 QUILL_LOG_WARNING(
447 log,
"Terminating dimer due to lost mode (align {:.3f}).", alignment);
449
450 if (bestNegativeCurvature < 0.0) {
451
452 C_tau = bestNegativeCurvature;
454 x0->setPositionsV(bestX0Positions);
455 x1->setPositionsV(bestX0Positions + delta * bestTau);
457 QUILL_LOG_DEBUG(
458 log,
"Restored best negative curvature state: C_tau={:.4f}",
C_tau);
459 throw eonc::DimerModeRestoredException();
460 } else {
461 QUILL_LOG_WARNING(
462 log,
"Never found negative curvature. Final C_tau: {:.4f}",
C_tau);
463 throw eonc::DimerModeLostException();
464 }
465 }
466
467 } while (std::abs(phi_prime) > std::abs(phi_tol) &&
468 std::abs(phi_min) > std::abs(phi_tol) &&
470}
bool foundNegativeCurvature
std::vector< VectorXd > s
VectorXd fixedReferenceMode
std::vector< VectorXd > positions
std::vector< VectorXd > gradients
std::vector< VectorXd > y
std::vector< double > rho
void rotationRemove(const AtomMatrix r1, std::shared_ptr< Matter > m2)
constexpr double safe_div(double num, double denom, double fallback=0.0)
double safe_acos(double x)
constexpr double safe_recip(double x, double fallback=0.0)
double safe_atan_ratio(double num, double denom, double fallback=0.0)
void evaluateTogether(Potential &pot, std::span< Matter *const > systems)
Evaluates every system that needs a force update.
bool potAllowsSharedInstance(const P &p) noexcept
std::optional< DimerRotationResult > runAlternativeRotation(DimerRotationBackend backend, const std::shared_ptr< Matter > &matter, const Parameters ¶ms, const std::shared_ptr< Potential > &pot, const AtomMatrix &initialDirection, quill::Logger *log=nullptr)
Run Lanczos, Davidson, or LOR.