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