Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2026 CP2K developers group <https://cp2k.org> */
4 : /* */
5 : /* SPDX-License-Identifier: GPL-2.0-or-later */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : #if defined(__ACE)
9 :
10 : // tested with lammps-user-pace-v.2023.11.25.fix2
11 :
12 : #if 0
13 : #include <stdio.h>
14 : #endif
15 : #include <string>
16 :
17 : #include "ace-evaluator/ace_c_basis.h"
18 : #include "ace-evaluator/ace_evaluator.h"
19 : #include "ace-evaluator/ace_recursive.h"
20 : #include "ace-evaluator/ace_version.h"
21 : #include "ace/ace_b_basis.h"
22 :
23 : struct ACEData {
24 : ACECTildeBasisSet *basis_set;
25 : ACERecursiveEvaluator *ace;
26 : Array1D<DOUBLE_TYPE> *virial;
27 : Array2D<DOUBLE_TYPE> *forces;
28 : };
29 :
30 6 : bool hasEnding(std::string const &fullString, std::string const &ending) {
31 6 : if (fullString.length() >= ending.length()) {
32 6 : return (0 == fullString.compare(fullString.length() - ending.length(),
33 6 : ending.length(), ending));
34 : } else {
35 : return false;
36 : }
37 : }
38 :
39 6 : extern "C" void AcePotInitialize(int ntypec, const char *symbolsc, int nlen,
40 : const char *potential_file_name, double *rcutc,
41 : void **acedata_ptr) {
42 :
43 : // avoid mixing C++ I/O with Fortran I/O, TODO: return this data so it can
44 : // be printed on the Fortran side.
45 : #if 0
46 : printf("---PACE initialization----\n");
47 :
48 : printf("ACE version: %d.%d.%d\n", VERSION_YEAR, VERSION_MONTH, VERSION_DAY);
49 : #endif
50 :
51 6 : std::string potential_file_name_str(potential_file_name);
52 : // trim potential_file_name
53 6 : potential_file_name_str.erase(
54 6 : potential_file_name_str.find_last_not_of(" \n\r\t") + 1);
55 :
56 6 : std::string symbols_str(symbolsc, ntypec * 2);
57 6 : std::vector<std::string> elements;
58 18 : for (int i = 0; i < ntypec; i++) {
59 12 : auto el_str = symbols_str.substr(2 * i, 2);
60 12 : el_str.erase(el_str.find_last_not_of(" \n\r\t") + 1);
61 12 : elements.push_back(el_str);
62 12 : }
63 6 : if (ntypec != elements.size())
64 0 : throw std::runtime_error(
65 0 : "Number of elements and elements list are inconsistent");
66 :
67 : #if 0
68 : printf("Number of atom types: %d\n", ntypec);
69 : printf("Element mapping: ");
70 : for (int i = 0; i < ntypec; i++)
71 : printf(" `%s`", elements[i].c_str());
72 : printf("\n");
73 : #endif
74 :
75 : // Elements are contained in a string of length 2*ntypec
76 : // Each element has two chars in the string
77 : // The sequence of the elements in the string corresponds to their mapping,
78 : // i.e., the first element in the string is element 1, the second element in
79 : // the string is element 2, the third element in the string is element 3,
80 : // ...
81 :
82 : #if 0
83 : printf("Filename: '%s'\n",
84 : potential_file_name_str.c_str());
85 : #endif
86 :
87 6 : ACEData *aceData;
88 :
89 6 : if (hasEnding(potential_file_name_str, ".yaml")) {
90 12 : ACEBBasisSet bBasisSet = ACEBBasisSet(potential_file_name_str);
91 6 : ACECTildeBasisSet cTildeBasisSet = bBasisSet.to_ACECTildeBasisSet();
92 6 : aceData = new ACEData;
93 6 : aceData->basis_set = new ACECTildeBasisSet(cTildeBasisSet);
94 6 : } else if (hasEnding(potential_file_name_str, ".yace")) {
95 0 : aceData = new ACEData;
96 0 : aceData->basis_set = new ACECTildeBasisSet(potential_file_name_str.c_str());
97 : } else {
98 0 : throw std::invalid_argument("Unrecognized file format: '" +
99 0 : potential_file_name_str + "'");
100 : }
101 6 : aceData->ace = new ACERecursiveEvaluator();
102 6 : aceData->ace->set_recursive(true);
103 :
104 6 : aceData->ace->element_type_mapping.init(1 + ntypec);
105 : // ace->element_type_mapping = {0,1,0}; // 0->0, 1(CP2k)-> 1(ACE),
106 : // 2(CP2k)-> 0(ACE)
107 18 : for (int i = 1; i <= ntypec; i++) {
108 12 : auto elemname = elements.at(i - 1);
109 12 : SPECIES_TYPE mu = aceData->basis_set->get_species_index_by_name(elemname);
110 12 : if (mu != -1) {
111 : #if 0
112 : printf("Mapping CP2K atom type #%d(%s) -> ACE species type #%d\n", i,
113 : elemname.c_str(), mu);
114 : #endif
115 : // set up CP2K atom type to ACE species mapping for ace evaluator
116 12 : aceData->ace->element_type_mapping(i) = mu;
117 : } else {
118 0 : throw std::runtime_error("Element " + elemname +
119 0 : " is not supported by ACE-potential from file " +
120 0 : potential_file_name_str);
121 : }
122 12 : }
123 :
124 : // the cutoffs of all pairs are stored in a ntypec*ntypec array
125 : // the index 0 corresponds to element 1
126 : // the index 1 corresponds to element 2
127 : // the index 2 corresponds to element 3, ...
128 : int k = 0;
129 18 : for (int i = 1; i <= ntypec; i++) {
130 36 : for (int j = 1; j <= ntypec; j++) {
131 24 : rcutc[k] = aceData->basis_set->radial_functions->cut(
132 24 : aceData->ace->element_type_mapping(i),
133 24 : aceData->ace->element_type_mapping(j));
134 24 : k++;
135 : }
136 : }
137 :
138 6 : aceData->ace->set_basis(*aceData->basis_set, 1);
139 6 : aceData->virial = new Array1D<DOUBLE_TYPE>(6, "virial");
140 6 : aceData->forces = new Array2D<DOUBLE_TYPE>(1, 3, "forces");
141 6 : *acedata_ptr = (void *)aceData;
142 : #if 0
143 : printf("---Done PACE initialization----\n");
144 : #endif
145 6 : }
146 :
147 6 : extern "C" void AcePotFinalize(void **acedata_ptr) {
148 6 : ACEData *aceData = (ACEData *)*acedata_ptr;
149 6 : delete aceData->basis_set;
150 6 : delete aceData->ace;
151 12 : delete aceData->virial;
152 12 : delete aceData->forces;
153 6 : delete aceData;
154 6 : }
155 :
156 206 : extern "C" void AcePotCompute(int natomc, int nghostc, int neic, int *neiatc,
157 : int *originc, int *nlistc, int *attypec,
158 : double *atposc, double *forcec, double *virialc,
159 : double *energyc, void **acedata_ptr) {
160 206 : ACEData *aceData = (ACEData *)*acedata_ptr;
161 : // re-point double **x (LAMMPS/C style of 2D array) to atposc
162 206 : int tot_nat = natomc + nghostc;
163 206 : double **x = new double *[tot_nat];
164 291357 : for (int i = 0; i < tot_nat; i++) {
165 291151 : x[i] = &atposc[3 * i];
166 : }
167 :
168 412 : std::vector<int> numneigh(natomc, 0);
169 39596 : for (int i = 0; i < natomc; i++) {
170 39390 : numneigh[i] = neiatc[i + 1] - neiatc[i];
171 : }
172 :
173 : // determine the maximum number of neighbours
174 : int i, jnum;
175 : int max_jnum = 0;
176 : int nei = 0;
177 39596 : for (i = 0; i < natomc; i++) {
178 39390 : jnum = numneigh[i];
179 39390 : nei = nei + jnum;
180 39390 : if (jnum > max_jnum)
181 : max_jnum = jnum;
182 : }
183 :
184 206 : aceData->ace->resize_neighbours_cache(max_jnum);
185 :
186 206 : double dx, dy, dz, fx, fy, fz;
187 :
188 : // resize forces array
189 206 : if (aceData->forces->get_dim(0) < natomc)
190 6 : aceData->forces->resize(natomc, 3);
191 :
192 206 : aceData->forces->fill(0);
193 206 : aceData->virial->fill(0);
194 :
195 : // main loop over atoms
196 39596 : for (i = 0; i < natomc; i++) {
197 39390 : jnum = numneigh[i];
198 :
199 39390 : const int *jlist = &nlistc[neiatc[i]];
200 :
201 39390 : aceData->ace->compute_atom(i, x, attypec, jnum, jlist);
202 :
203 39390 : energyc[i] = aceData->ace->e_atom;
204 :
205 39390 : const double xtmp = x[i][0];
206 39390 : const double ytmp = x[i][1];
207 39390 : const double ztmp = x[i][2];
208 :
209 4260152 : for (int jj = 0; jj < jnum; jj++) {
210 4220762 : int j = jlist[jj];
211 :
212 4220762 : dx = x[j][0] - xtmp;
213 4220762 : dy = x[j][1] - ytmp;
214 4220762 : dz = x[j][2] - ztmp;
215 :
216 4220762 : fx = aceData->ace->neighbours_forces(jj, 0);
217 4220762 : fy = aceData->ace->neighbours_forces(jj, 1);
218 4220762 : fz = aceData->ace->neighbours_forces(jj, 2);
219 :
220 4220762 : (*aceData->forces)(i, 0) += fx;
221 4220762 : (*aceData->forces)(i, 1) += fy;
222 4220762 : (*aceData->forces)(i, 2) += fz;
223 :
224 : // virial f_dot_r, identical to LAMMPS virial_fdotr_compute
225 4220762 : (*aceData->virial)(0) += dx * fx;
226 4220762 : (*aceData->virial)(1) += dy * fy;
227 4220762 : (*aceData->virial)(2) += dz * fz;
228 4220762 : (*aceData->virial)(3) += dx * fy;
229 4220762 : (*aceData->virial)(4) += dx * fz;
230 4220762 : (*aceData->virial)(5) += dy * fz;
231 :
232 : // update forces only for real atoms
233 4220762 : if (j < natomc) {
234 1723464 : (*aceData->forces)(j, 0) -= fx;
235 1723464 : (*aceData->forces)(j, 1) -= fy;
236 1723464 : (*aceData->forces)(j, 2) -= fz;
237 : } else {
238 : // map ghost j into true_j within periodic cell
239 2497298 : int true_j = originc[j];
240 2497298 : (*aceData->forces)(true_j, 0) -= fx;
241 2497298 : (*aceData->forces)(true_j, 1) -= fy;
242 2497298 : (*aceData->forces)(true_j, 2) -= fz;
243 : }
244 : }
245 : }
246 :
247 : double ene = 0.0;
248 39596 : for (int i = 0; i < natomc; i++) {
249 39390 : ene += energyc[i];
250 :
251 39390 : forcec[3 * i] = (*aceData->forces)(i, 0);
252 39390 : forcec[3 * i + 1] = (*aceData->forces)(i, 1);
253 39390 : forcec[3 * i + 2] = (*aceData->forces)(i, 2);
254 : }
255 :
256 : // copy virials
257 1442 : for (int i = 0; i < 6; i++) {
258 1236 : virialc[i] = (*aceData->virial)(i);
259 : }
260 :
261 206 : delete[] x;
262 206 : }
263 :
264 : #endif // defined(__ACE)
265 :
266 : // EOF
|
