Loading...
Searching...
No Matches
eonc::ARTnSaddleSearch Class Reference

Saddle search method using the Activation-Relaxation Technique nouveau. More...

#include <ARTnSaddleSearch.h>

Inheritance diagram for eonc::ARTnSaddleSearch:

Public Member Functions

 ARTnSaddleSearch (std::shared_ptr< Matter > matterPassed, std::shared_ptr< Potential > potPassed, AtomMatrix modeInitial, const Parameters &paramsPassed)
 ~ARTnSaddleSearch () override
int run () override
double getEigenvalue () override
AtomMatrix getEigenvector () override
std::string_view describeStatus (int status) const override
int getStatus () const override
int getIterationCount () const override
int getForceCalls () const override
Public Member Functions inherited from eonc::SaddleSearchMethod
 SaddleSearchMethod (std::shared_ptr< Potential > potPassed, const Parameters &paramsPassed)
virtual ~SaddleSearchMethod ()

Static Public Attributes

static constexpr int STATUS_GOOD = 0
static constexpr int STATUS_BAD_MAX_ITERATIONS
static constexpr int STATUS_BAD_ARTN_ERROR = 22

Private Attributes

std::shared_ptr< Mattermatter
double eigenvalue {std::numeric_limits<double>::quiet_NaN()}
AtomMatrix eigenvector
AtomMatrix mode
int status {0}
int iteration {0}
int forcecalls {0}
eonc::log::Scoped log

Additional Inherited Members

Protected Attributes inherited from eonc::SaddleSearchMethod
std::shared_ptr< Potentialpot
const Parametersparams

Detailed Description

Saddle search method using the Activation-Relaxation Technique nouveau.

Wraps the pARTn Fortran library via its C API (artn.h).

Definition at line 47 of file ARTnSaddleSearch.h.

Constructor & Destructor Documentation

◆ ARTnSaddleSearch()

eonc::ARTnSaddleSearch::ARTnSaddleSearch ( std::shared_ptr< Matter > matterPassed,
std::shared_ptr< Potential > potPassed,
AtomMatrix modeInitial,
const Parameters & paramsPassed )

Definition at line 21 of file ARTnSaddleSearch.cpp.

25 : SaddleSearchMethod(potPassed, paramsPassed),
26 matter{matterPassed},
27 mode{modeInitial},
28 eigenvector{AtomMatrix::Zero(matterPassed->numberOfAtoms(), 3)} {
30 if (!log) {
31 throw std::runtime_error("ARTnSaddleSearch: Logger not initialized");
32 }
33}
eonc::log::Scoped log
std::shared_ptr< Matter > matter
SaddleSearchMethod(std::shared_ptr< Potential > potPassed, const Parameters &paramsPassed)
quill::Logger * get() noexcept
Get or create the default "combi" logger.
Definition EonLogger.h:44

◆ ~ARTnSaddleSearch()

eonc::ARTnSaddleSearch::~ARTnSaddleSearch ( )
override

Definition at line 35 of file ARTnSaddleSearch.cpp.

35 {
36#ifdef WITH_ARTN
37 // Clean up is done within the search loop, not in destructor
38#endif
39}

Member Function Documentation

◆ describeStatus()

std::string_view eonc::ARTnSaddleSearch::describeStatus ( int status) const
overridevirtual

Implements eonc::SaddleSearchMethod.

Definition at line 441 of file ARTnSaddleSearch.cpp.

441 {
442 switch (status) {
443 case STATUS_GOOD:
444 return "Success";
446 return "Too many iterations";
448 return "ARTn backend error";
449 default:
450 return "Unknown status";
451 }
452}
static constexpr int STATUS_GOOD
static constexpr int STATUS_BAD_ARTN_ERROR
static constexpr int STATUS_BAD_MAX_ITERATIONS

◆ getEigenvalue()

double eonc::ARTnSaddleSearch::getEigenvalue ( )
overridevirtual

Implements eonc::SaddleSearchMethod.

Definition at line 427 of file ARTnSaddleSearch.cpp.

427 {
428 if (log && std::isnan(eigenvalue)) {
429 QUILL_LOG_WARNING(log, "Requesting uninitialized/invalid eigenvalue");
430 }
431 return eigenvalue;
432}

◆ getEigenvector()

AtomMatrix eonc::ARTnSaddleSearch::getEigenvector ( )
overridevirtual

Implements eonc::SaddleSearchMethod.

Definition at line 434 of file ARTnSaddleSearch.cpp.

434 {
435 if (log && std::isnan(eigenvalue)) {
436 QUILL_LOG_WARNING(log, "Requesting uninitialized eigenvector");
437 }
438 return eigenvector;
439}

◆ getForceCalls()

int eonc::ARTnSaddleSearch::getForceCalls ( ) const
inlineoverridevirtual

Reimplemented from eonc::SaddleSearchMethod.

Definition at line 65 of file ARTnSaddleSearch.h.

◆ getIterationCount()

int eonc::ARTnSaddleSearch::getIterationCount ( ) const
inlineoverridevirtual

Reimplemented from eonc::SaddleSearchMethod.

Definition at line 64 of file ARTnSaddleSearch.h.

◆ getStatus()

int eonc::ARTnSaddleSearch::getStatus ( ) const
inlineoverridevirtual

Reimplemented from eonc::SaddleSearchMethod.

Definition at line 63 of file ARTnSaddleSearch.h.

63{ return status; }

◆ run()

int eonc::ARTnSaddleSearch::run ( void )
overridevirtual

Implements eonc::SaddleSearchMethod.

Definition at line 41 of file ARTnSaddleSearch.cpp.

41 {
42#ifdef WITH_ARTN
43 auto &res = get_artn_resource();
44 const int nat = matter->numberOfAtoms();
45
46 if (mode.rows() != nat || mode.cols() != 3) {
47 mode = AtomMatrix::Zero(nat, 3);
48 }
49
50 // Pre-declare force-loop storage outside the lock so it outlives the
51 // setup critical section. These reads do not touch pARTn state.
52 AtomMatrix positions = matter->getPositions();
53 AtomMatrix forces = AtomMatrix::Zero(nat, 3);
54 AtomMatrix displacement = AtomMatrix::Zero(nat, 3);
55
56 // Eigen::Maps give Fortran a column-major [3, nat] view over the same
57 // memory as the row-major [nat, 3] AtomMatrix (zero-copy).
58 Eigen::Map<AtomMatrixF> pos_map(positions.data(), 3, nat);
59 Eigen::Map<AtomMatrixF> force_map(forces.data(), 3, nat);
60 Eigen::Map<AtomMatrixF> disp_map(displacement.data(), 3, nat);
61 Eigen::Map<AtomMatrixF> mode_map(mode.data(), 3, nat);
62
63 const double push_step = params.artn_options.push_step_size;
64 const double mode_norm = mode.norm();
65 AtomMatrixF mode_fort;
66 if (mode_norm > 1e-10) {
67 mode_fort = mode_map;
68 Eigen::Map<VectorXd> mode_vec_map(mode_fort.data(), mode_fort.size());
69 mode_vec_map *= (push_step / mode_norm);
70 }
71 int dim_mode[2] = {3, nat};
72
73 // 1. Library Initialization, Configuration, and Initial Push (Locked)
74 // Held as a single critical section so concurrent ARTn searches in
75 // the same process cannot interleave create(), set_param(), setup,
76 // or push_init calls on pARTn's non-thread-safe global state.
77 {
78 std::lock_guard<std::mutex> lock(res.library_mutex);
79
80 try {
81 res.require_loaded();
82 } catch (const std::exception &e) {
83 QUILL_LOG_ERROR(log, "ARTn library not available: {}", e.what());
85 return status;
86 }
87
88 res.get_create_fn()();
89
90 // Set engine units first (required before any other params)
91 const char *units = "lammps/metal";
92 int size0 = 0;
93 int result_units = res.get_set_param_fn()("engine_units", 0, &size0, units);
94 if (result_units != 0) {
95 QUILL_LOG_ERROR(log, "set_param(engine_units) failed with code {}",
96 result_units);
97 }
98
99 // Set parameters from eOn config
100 double push_step = params.artn_options.push_step_size;
101 int result_push =
102 res.get_set_param_fn()("push_step_size", 0, &size0, &push_step);
103 if (result_push != 0) {
104 QUILL_LOG_ERROR(log, "set_param(push_step_size) failed with code {}",
105 result_push);
106 }
107
108 double force_thr = params.artn_options.force_threshold;
109 int result_force =
110 res.get_set_param_fn()("forc_thr", 0, &size0, &force_thr);
111 if (result_force != 0) {
112 QUILL_LOG_ERROR(log, "set_param(forc_thr) failed with code {}",
113 result_force);
114 }
115
116 // filin names the artn.in input file pARTn reads at setup. artn_create
117 // resets its internal value to NAN_STR ("BBBB") meaning "undefined", so
118 // an empty eOn config leaves pARTn reading no file at all. If the user
119 // does set a path, surface a missing file before setup_artn runs so the
120 // failure names the file instead of hiding inside pARTn's ERR_FILE code.
121 const std::string &filin = params.artn_options.filin;
122 if (!filin.empty()) {
123 if (!std::filesystem::exists(filin)) {
124 QUILL_LOG_ERROR(log, "artn_options.filin '{}' does not exist", filin);
125 res.get_destroy_fn()();
127 return status;
128 }
129 int result_filin =
130 res.get_set_param_fn()("filin", 0, &size0, filin.c_str());
131 if (result_filin != 0) {
132 QUILL_LOG_ERROR(log, "set_param(filin) failed with code {}",
133 result_filin);
134 }
135 }
136
137 int verbosity = 3;
138 int result_verbose =
139 res.get_set_param_fn()("verbose", 0, &size0, &verbosity);
140 if (result_verbose != 0) {
141 QUILL_LOG_ERROR(log, "set_param(verbose) failed with code {}",
142 result_verbose);
143 }
144
145 // ninit controls initial push steps before Lanczos eigenmode estimation.
146 // 0 = skip push, go straight to Lanczos (appropriate when eOn provides
147 // the displacement direction via push_init); >0 = push that many steps.
148 // -1 sentinel means "leave pARTn's own default in place", so we only
149 // call set_param when the user asked for a specific value.
150 if (params.artn_options.ninit >= 0) {
151 int ninit = params.artn_options.ninit;
152 int result_ninit = res.get_set_param_fn()("ninit", 0, &size0, &ninit);
153 if (result_ninit != 0) {
154 QUILL_LOG_ERROR(log, "set_param(ninit) failed with code {}",
155 result_ninit);
156 }
157 }
158
159 // nperp_limitation: controls perp-relax steps per Lanczos cycle.
160 // pARTn defaults are tuned for exploration from minimum. For refinement
161 // near a saddle, -1 (unlimited) or 20-30 (for ML potentials) is better.
162 if (params.artn_options.nperp_limitation != "default") {
163 // Parse comma-separated integers into a vector
164 std::vector<int> nperp_vals;
165 std::istringstream ss(params.artn_options.nperp_limitation);
166 std::string token;
167 while (std::getline(ss, token, ',')) {
168 nperp_vals.push_back(std::stoi(token));
169 }
170 if (!nperp_vals.empty()) {
171 int nperp_size = static_cast<int>(nperp_vals.size());
172 int result_nperp = res.get_set_param_fn()(
173 "nperp_limitation", 1, &nperp_size, nperp_vals.data());
174 if (result_nperp != 0) {
175 QUILL_LOG_WARNING(log, "set_param(nperp_limitation) failed: {}",
176 result_nperp);
177 }
178 }
179 }
180
181 // lanczos_min_size: minimum Lanczos iterations before convergence check.
182 // Default 3 for exploration; 1 for refinement near a saddle.
183 if (params.artn_options.lanczos_min_size >= 0) {
184 int lms = params.artn_options.lanczos_min_size;
185 res.get_set_param_fn()("lanczos_min_size", 0, &size0, &lms);
186 }
187
188 // nsmooth: number of smooth interpolation steps. 0 disables.
189 if (params.artn_options.nsmooth >= 0) {
190 int ns = params.artn_options.nsmooth;
191 res.get_set_param_fn()("nsmooth", 0, &size0, &ns);
192 }
193
194 // nnewchance: retries permitted when Lanczos returns a positive lowest
195 // eigenvalue (convex region, no unstable mode). pARTn defaults to 0,
196 // i.e. immediate "EIGENVALUE LOST" failure -- too brittle for small
197 // clusters where the eigenvalue can flip positive transiently before
198 // the saddle direction stabilizes. eOn's default of 3 (Parameters.h)
199 // gives Lanczos a few random-restart chances before giving up.
200 if (params.artn_options.nnewchance >= 0) {
201 int nnc = params.artn_options.nnewchance;
202 res.get_set_param_fn()("nnewchance", 0, &size0, &nnc);
203 }
204
205 // Setup
206 bool cerr = false;
207 res.get_setup_fn()(nat, &cerr);
208 if (cerr) {
209 QUILL_LOG_ERROR(log, "ARTn setup failed (nat={})", nat);
210 res.get_destroy_fn()();
212 return status;
213 }
214
215 // Initial push vector (if eOn supplied a non-trivial mode). Must be set
216 // after setup_artn() and before the first artn_step(); keep inside the
217 // same critical section so no concurrent ARTn search can re-init between
218 // setup and push.
219 if (mode_norm > 1e-10) {
220 int result =
221 res.get_set_param_fn()("push_init", 2, dim_mode, mode_fort.data());
222 if (result != 0) {
223 QUILL_LOG_WARNING(log, "set_param(push_init) failed with code {}",
224 result);
225 }
226 }
227 }
228
229 // Per-atom metadata for the Fortran step (no pARTn state, unlocked).
230 std::vector<int> ityp(nat);
231 std::vector<int> if_pos(3 * nat, 1); // all atoms free by default
232 double box_f[9];
233 bool lconv = false;
234
235 for (int i = 0; i < nat; i++) {
236 if (matter->getFixed(i)) {
237 Eigen::Map<Eigen::Vector3i>(&if_pos[i * 3]).setZero();
238 }
239 ityp[i] = matter->getAtomicNr(i);
240 }
241
242 // Convert box to column-major 3x3 - copy to array
243 Matrix3d cell = matter->getCell();
244 for (int i = 0; i < 3; i++)
245 for (int j = 0; j < 3; j++)
246 box_f[j * 3 + i] = cell(i, j);
247
248 int maxIter = params.artn_options.max_iterations;
249
250 // while-loop (not for) so `iteration` reports the index of the converged
251 // step rather than one past it: on convergence during step k, we break
252 // before the increment and the post-loop value is k, matching log output.
253 iteration = 0;
254 while (iteration < maxIter && !lconv) {
255 // 2. Parallel PES Evaluation (UNLOCKED)
256 // Concurrent ARTn searches in the same process would share the same
257 // PES call, so the potential evaluation itself is not serialized here.
258 double energy = matter->getPotentialEnergy();
259 forces = matter->getForces(); // Returns RowMajor Nx3
260 this->forcecalls++;
261
262 // 3. ARTn State Update (LOCKED)
263 // Serialize only the interaction with the non-thread-safe backend.
264 //
265 // Perf note (artn-plugin >= 9dab2053): the inner Lanczos eigenvector
266 // reconstruction now uses intrinsic matmul on a reshaped Vmat slice,
267 // which allocates a temporary [3*nat, ilanc] array per Lanczos iteration.
268 // Negligible at our sizes (small molecules, ilanc < O(20)); revisit if
269 // we ever drive artn against large supercell DFT.
270 {
271 std::unique_lock<std::mutex> lock(res.library_mutex);
272 res.get_artn_step_fn()(nat, energy, force_map.data(), ityp.data(),
273 pos_map.data(), box_f, if_pos.data(),
274 disp_map.data(), &lconv);
275 }
276
277 if (!lconv) {
278 // Apply displacement and update Matter (PES-specific, typically
279 // thread-safe) Add displacement to positions
280 positions += displacement;
281 matter->setPositions(positions);
282 iteration++;
283 }
284 }
285
286 // 4. Data Retrieval (LOCKED)
287 if (lconv) {
288 std::lock_guard<std::mutex> lock(res.library_mutex);
289
290 // pARTn exposes a dedicated C get_error() that returns both the error
291 // code and a c_malloc'd message pointer (see m_artn_error.f90 in
292 // artn-plugin). Prefer it when available; fall back to the has_error
293 // flag via get_data for older libartn builds that predate the wrapper.
294 int artn_err = 0;
295 std::string artn_err_msg;
296 if (auto *get_error_fn_ = res.get_get_error_fn()) {
297 void *cmsg = nullptr;
298 artn_err = get_error_fn_(&cmsg);
299 if (artn_err != 0 && cmsg != nullptr) {
300 artn_err_msg.assign(static_cast<const char *>(cmsg));
301 std::free(cmsg);
302 }
303 } else {
304 bool *has_error_ptr = nullptr;
305 int result_has_error = res.get_get_data_fn()(
306 "has_error", reinterpret_cast<void **>(&has_error_ptr));
307 if (result_has_error != 0) {
308 QUILL_LOG_WARNING(log, "get_data(has_error) failed with code {}",
309 result_has_error);
310 } else if (has_error_ptr) {
311 artn_err = *has_error_ptr ? 1 : 0;
312 std::free(has_error_ptr);
313 }
314 }
315 bool has_error = artn_err != 0;
316
317 bool has_sad = false;
318 bool *has_sad_ptr = nullptr;
319 int result_has_sad = res.get_get_data_fn()(
320 "has_sad", reinterpret_cast<void **>(&has_sad_ptr));
321 if (result_has_sad != 0) {
322 QUILL_LOG_WARNING(log, "get_data(has_sad) failed with code {}",
323 result_has_sad);
324 } else if (has_sad_ptr) {
325 has_sad = *has_sad_ptr;
326 std::free(has_sad_ptr);
327 }
328
329 if (has_error && !artn_err_msg.empty()) {
330 QUILL_LOG_WARNING(log, "pARTn reported error {}: {}", artn_err,
331 artn_err_msg);
332 }
333
334 // If a saddle was found, accept it even if force didn't fully converge
335 if (has_sad) {
336 QUILL_LOG_INFO(
337 log,
338 "ARTn found saddle after {} iterations (has_error={}, has_sad={})",
339 this->iteration, has_error, has_sad);
340
341 // Retrieve the saddle coordinates tracked internally by pARTn so the
342 // Matter object matches the reported eigenpair and subsequent endpoint
343 // minimizations start from the actual saddle.
344 double *tau_sad_ptr = nullptr;
345 int result_tau_sad = res.get_get_data_fn()(
346 "tau_sad", reinterpret_cast<void **>(&tau_sad_ptr));
347 if (result_tau_sad == 0 && tau_sad_ptr) {
349 std::vector<double>(tau_sad_ptr, tau_sad_ptr + 3 * nat), nat));
350 std::free(tau_sad_ptr);
351 } else {
352 QUILL_LOG_WARNING(
353 log, "Failed to retrieve tau_sad (result={}, ptr_valid={})",
354 result_tau_sad, tau_sad_ptr != nullptr);
355 }
356
357 // Retrieve eigenvalue
358 double *eigval_ptr = nullptr;
359 int result_eigval = res.get_get_data_fn()(
360 "eigval_sad", reinterpret_cast<void **>(&eigval_ptr));
361 if (result_eigval == 0 && eigval_ptr) {
362 eigenvalue = *eigval_ptr;
363 std::free(eigval_ptr);
364 } else {
365 QUILL_LOG_WARNING(
366 log, "Failed to retrieve eigenvalue (result={}, ptr_valid={})",
367 result_eigval, eigval_ptr != nullptr);
368 eigenvalue =
369 std::numeric_limits<double>::quiet_NaN(); // Use NaN to indicate
370 // missing value
371 }
372
373 // Retrieve eigenvector (3*nat flat array, column-major from Fortran).
374 // get_data allocates via c_malloc and writes the pointer to cval.
375 // The C header says void* but the Fortran intent(out) semantics
376 // require void** (see artn_c_wrappers.f90:324 and LAMMPS example).
377 double *evec_ptr = nullptr;
378 int result_evec = res.get_get_data_fn()(
379 "eigen_sad", reinterpret_cast<void **>(&evec_ptr));
380 if (result_evec == 0 && evec_ptr) {
381 // Use direct Eigen::Map to convert from Fortran layout
383 std::vector<double>(evec_ptr, evec_ptr + 3 * nat), nat);
384
385 // get_data allocates via c_malloc (artn_c_wrappers.f90), safe to free
386 std::free(evec_ptr);
387 } else {
388 QUILL_LOG_ERROR(
389 log, "Failed to retrieve eigenvector (result={}, ptr_valid={})",
390 result_evec, evec_ptr != nullptr);
391 // Set eigenvector to zero matrix if retrieval failed
392 eigenvector = AtomMatrix::Zero(nat, 3);
393 }
394
396 res.get_destroy_fn()();
397 return status;
398 }
399
400 // No saddle found - this is a real error
401 QUILL_LOG_WARNING(
402 log, "ARTn stopped after {} iterations (has_error={}, has_sad={})",
403 iteration, has_error, has_sad);
405 res.get_destroy_fn()();
406 return status;
407 }
408
409 QUILL_LOG_WARNING(log, "ARTn did not converge after {} iterations",
410 iteration);
412
413 // Clean up in all cases
414 {
415 std::lock_guard<std::mutex> lock(res.library_mutex);
416 res.get_destroy_fn()();
417 }
418 return status;
419
420#else
421 QUILL_LOG_ERROR(log, "ARTn support not compiled");
423 return status;
424#endif
425}
Eigen::Matrix< double, 3, 3, eOnStorageOrder > Matrix3d
Definition Eigen.h:35
Eigen::Matrix< double, 3, Eigen::Dynamic, Eigen::ColMajor > AtomMatrixF
Definition Eigen.h:44
Eigen::Matrix< double, Eigen::Dynamic, 3, eOnStorageOrder > AtomMatrix
Definition Eigen.h:37
ARTnResource & get_artn_resource()
Global access to thread-safe ARTn resource.
AtomMatrix from_fortran_layout_vector(const std::vector< double > &flat_colmajor, int nat)
Reconstruct AtomMatrix from a flat column-major vector (e.g.
Definition Eigen.h:61

Member Data Documentation

◆ eigenvalue

double eonc::ARTnSaddleSearch::eigenvalue {std::numeric_limits<double>::quiet_NaN()}
private

Definition at line 69 of file ARTnSaddleSearch.h.

69{std::numeric_limits<double>::quiet_NaN()};

◆ eigenvector

AtomMatrix eonc::ARTnSaddleSearch::eigenvector
private

Definition at line 70 of file ARTnSaddleSearch.h.

◆ forcecalls

int eonc::ARTnSaddleSearch::forcecalls {0}
private

Definition at line 73 of file ARTnSaddleSearch.h.

73{0};

◆ iteration

int eonc::ARTnSaddleSearch::iteration {0}
private

Definition at line 72 of file ARTnSaddleSearch.h.

72{0};

◆ log

eonc::log::Scoped eonc::ARTnSaddleSearch::log
private

Definition at line 74 of file ARTnSaddleSearch.h.

◆ matter

std::shared_ptr<Matter> eonc::ARTnSaddleSearch::matter
private

Definition at line 68 of file ARTnSaddleSearch.h.

◆ mode

AtomMatrix eonc::ARTnSaddleSearch::mode
private

Definition at line 70 of file ARTnSaddleSearch.h.

◆ status

int eonc::ARTnSaddleSearch::status {0}
private

Definition at line 71 of file ARTnSaddleSearch.h.

71{0};

◆ STATUS_BAD_ARTN_ERROR

int eonc::ARTnSaddleSearch::STATUS_BAD_ARTN_ERROR = 22
staticconstexpr

Definition at line 52 of file ARTnSaddleSearch.h.

◆ STATUS_BAD_MAX_ITERATIONS

int eonc::ARTnSaddleSearch::STATUS_BAD_MAX_ITERATIONS
staticconstexpr

◆ STATUS_GOOD

int eonc::ARTnSaddleSearch::STATUS_GOOD = 0
staticconstexpr

Definition at line 49 of file ARTnSaddleSearch.h.


The documentation for this class was generated from the following files: