41 {
42#ifdef WITH_ARTN
44 const int nat =
matter->numberOfAtoms();
45
46 if (
mode.rows() != nat ||
mode.cols() != 3) {
47 mode = AtomMatrix::Zero(nat, 3);
48 }
49
50
51
54 AtomMatrix displacement = AtomMatrix::Zero(nat, 3);
55
56
57
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();
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
74
75
76
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());
86 }
87
88 res.get_create_fn()();
89
90
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
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
117
118
119
120
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()();
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
146
147
148
149
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
160
161
162 if (
params.artn_options.nperp_limitation !=
"default") {
163
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
182
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
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
195
196
197
198
199
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
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()();
213 }
214
215
216
217
218
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
230 std::vector<int> ityp(nat);
231 std::vector<int> if_pos(3 * nat, 1);
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
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
251
252
255
256
257
258 double energy =
matter->getPotentialEnergy();
259 forces =
matter->getForces();
261
262
263
264
265
266
267
268
269
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
279
280 positions += displacement;
281 matter->setPositions(positions);
283 }
284 }
285
286
287 if (lconv) {
288 std::lock_guard<std::mutex> lock(res.library_mutex);
289
290
291
292
293
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
335 if (has_sad) {
336 QUILL_LOG_INFO(
338 "ARTn found saddle after {} iterations (has_error={}, has_sad={})",
340
341
342
343
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
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) {
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);
369 std::numeric_limits<double>::quiet_NaN();
370
371 }
372
373
374
375
376
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
383 std::vector<double>(evec_ptr, evec_ptr + 3 * nat), nat);
384
385
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
393 }
394
396 res.get_destroy_fn()();
398 }
399
400
401 QUILL_LOG_WARNING(
402 log,
"ARTn stopped after {} iterations (has_error={}, has_sad={})",
405 res.get_destroy_fn()();
407 }
408
409 QUILL_LOG_WARNING(
log,
"ARTn did not converge after {} iterations",
412
413
414 {
415 std::lock_guard<std::mutex> lock(res.library_mutex);
416 res.get_destroy_fn()();
417 }
419
420#else
421 QUILL_LOG_ERROR(
log,
"ARTn support not compiled");
424#endif
425}
Eigen::Matrix< double, 3, 3, eOnStorageOrder > Matrix3d
Eigen::Matrix< double, 3, Eigen::Dynamic, Eigen::ColMajor > AtomMatrixF
Eigen::Matrix< double, Eigen::Dynamic, 3, eOnStorageOrder > AtomMatrix
const Parameters & params
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.