1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 /* ----------------------------------------------------------------------
16 Contributing authors: Axel Kohlmeyer (Temple U), Venkatesh Botu, James Chapman (Lawrence Livermore National Lab)
17 ------------------------------------------------------------------------- */
18
19 #include "pair_agni.h"
20
21 #include "atom.h"
22 #include "citeme.h"
23 #include "comm.h"
24 #include "force.h"
25 #include "error.h"
26 #include "math_const.h"
27 #include "math_special.h"
28 #include "memory.h"
29 #include "neigh_list.h"
30 #include "neigh_request.h"
31 #include "neighbor.h"
32 #include "tokenizer.h"
33 #include "potential_file_reader.h"
34
35 #include <cmath>
36 #include <cstring>
37
38 using namespace LAMMPS_NS;
39 using namespace MathSpecial;
40
41 static const char cite_pair_agni[] =
42 "pair agni command:\n\n"
43 "@article{huan2019jpc,\n"
44 " author = {Huan, T. and Batra, R. and Chapman, J. and Kim, C. and Chandrasekaran, A. and Ramprasad, Rampi},\n"
45 " journal = {J. Phys. Chem. C},\n"
46 " volume = {121},\n"
47 " number = {34},\n"
48 " pages = {20715},\n"
49 " year = {2019},\n"
50 "}\n\n";
51
52 #define MAXLINE 10240
53 #define MAXWORD 40
54
55 /* ---------------------------------------------------------------------- */
56
PairAGNI(LAMMPS * lmp)57 PairAGNI::PairAGNI(LAMMPS *lmp) : Pair(lmp)
58 {
59 if (lmp->citeme) lmp->citeme->add(cite_pair_agni);
60
61 single_enable = 0;
62 restartinfo = 0;
63 one_coeff = 1;
64 manybody_flag = 1;
65 atomic_feature_version = 0;
66
67 centroidstressflag = CENTROID_NOTAVAIL;
68
69 no_virial_fdotr_compute = 1;
70
71 params = nullptr;
72 cutmax = 0.0;
73 }
74
75 /* ----------------------------------------------------------------------
76 check if allocated, since class can be destructed when incomplete
77 ------------------------------------------------------------------------- */
78
~PairAGNI()79 PairAGNI::~PairAGNI()
80 {
81 if (params) {
82 for (int i = 0; i < nparams; ++i) {
83 memory->destroy(params[i].eta);
84 memory->destroy(params[i].alpha);
85 memory->destroy(params[i].xU);
86 }
87 memory->destroy(params);
88 params = nullptr;
89 }
90 memory->destroy(elem1param);
91
92 if (allocated) {
93 memory->destroy(setflag);
94 memory->destroy(cutsq);
95 }
96 }
97
98 /* ---------------------------------------------------------------------- */
99
compute(int eflag,int vflag)100 void PairAGNI::compute(int eflag, int vflag)
101 {
102 int i,j,k,ii,jj,inum,jnum,itype;
103 double xtmp,ytmp,ztmp,delx,dely,delz;
104 double rsq;
105 int *ilist,*jlist,*numneigh,**firstneigh;
106
107 ev_init(eflag,vflag);
108
109 double **x = atom->x;
110 double **f = atom->f;
111 int *type = atom->type;
112
113 inum = list->inum;
114 ilist = list->ilist;
115 numneigh = list->numneigh;
116 firstneigh = list->firstneigh;
117
118 double fxtmp,fytmp,fztmp;
119 double *Vx, *Vy, *Vz;
120
121 // loop over full neighbor list of my atoms
122 for (ii = 0; ii < inum; ii++) {
123 i = ilist[ii];
124 itype = map[type[i]];
125 xtmp = x[i][0];
126 ytmp = x[i][1];
127 ztmp = x[i][2];
128 fxtmp = fytmp = fztmp = 0.0;
129
130 const Param &iparam = params[elem1param[itype]];
131 Vx = new double[iparam.numeta];
132 Vy = new double[iparam.numeta];
133 Vz = new double[iparam.numeta];
134 memset(Vx,0,iparam.numeta*sizeof(double));
135 memset(Vy,0,iparam.numeta*sizeof(double));
136 memset(Vz,0,iparam.numeta*sizeof(double));
137
138 jlist = firstneigh[i];
139 jnum = numneigh[i];
140
141 for (jj = 0; jj < jnum; jj++) {
142 j = jlist[jj];
143 j &= NEIGHMASK;
144
145 delx = xtmp - x[j][0];
146 dely = ytmp - x[j][1];
147 delz = ztmp - x[j][2];
148 rsq = delx*delx + dely*dely + delz*delz;
149
150 if ((rsq > 0.0) && (rsq < iparam.cutsq)) {
151 const double r = sqrt(rsq);
152 const double cF = 0.5*(cos((MathConst::MY_PI*r)/iparam.cut)+1.0);
153 const double wX = cF*delx/r;
154 const double wY = cF*dely/r;
155 const double wZ = cF*delz/r;
156
157 for (k = 0; k < iparam.numeta; ++k) {
158 double e = 0.0;
159
160 if (atomic_feature_version == AGNI_VERSION_1)
161 e = exp(-(iparam.eta[k]*rsq));
162 else if (atomic_feature_version == AGNI_VERSION_2)
163 e = (1.0 / (square(iparam.eta[k]) * iparam.gwidth * sqrt(MathConst::MY_2PI))) * exp(-(square(r - iparam.eta[k])) / (2.0 * square(iparam.gwidth)));
164
165 Vx[k] += wX*e;
166 Vy[k] += wY*e;
167 Vz[k] += wZ*e;
168 }
169 }
170 }
171
172 for (j = 0; j < iparam.numtrain; ++j) {
173 double kx = 0.0;
174 double ky = 0.0;
175 double kz = 0.0;
176
177 for (int k = 0; k < iparam.numeta; ++k) {
178 const double xu = iparam.xU[k][j];
179
180 kx += square(Vx[k] - xu);
181 ky += square(Vy[k] - xu);
182 kz += square(Vz[k] - xu);
183 }
184 const double e = -0.5/(square(iparam.sigma));
185 fxtmp += iparam.alpha[j]*exp(kx*e);
186 fytmp += iparam.alpha[j]*exp(ky*e);
187 fztmp += iparam.alpha[j]*exp(kz*e);
188 }
189 fxtmp += iparam.b;
190 fytmp += iparam.b;
191 fztmp += iparam.b;
192 f[i][0] += fxtmp;
193 f[i][1] += fytmp;
194 f[i][2] += fztmp;
195
196 if (evflag) ev_tally_xyz_full(i,0.0,0.0,fxtmp,fytmp,fztmp,delx,dely,delz);
197
198 delete [] Vx;
199 delete [] Vy;
200 delete [] Vz;
201 }
202
203 if (vflag_fdotr) virial_fdotr_compute();
204 }
205
206 /* ---------------------------------------------------------------------- */
207
allocate()208 void PairAGNI::allocate()
209 {
210 allocated = 1;
211 int n = atom->ntypes;
212
213 memory->create(setflag,n+1,n+1,"pair:setflag");
214 memory->create(cutsq,n+1,n+1,"pair:cutsq");
215 map = new int[n+1];
216 }
217
218 /* ----------------------------------------------------------------------
219 global settings
220 ------------------------------------------------------------------------- */
221
settings(int narg,char **)222 void PairAGNI::settings(int narg, char ** /* arg */)
223 {
224 if (narg != 0) error->all(FLERR,"Illegal pair_style command");
225 }
226
227 /* ----------------------------------------------------------------------
228 set coeffs for one or more type pairs
229 ------------------------------------------------------------------------- */
230
coeff(int narg,char ** arg)231 void PairAGNI::coeff(int narg, char **arg)
232 {
233 if (!allocated) allocate();
234
235 map_element2type(narg-3,arg+3);
236
237 if (nelements != 1)
238 error->all(FLERR,"Cannot handle multi-element systems with this potential");
239
240 // read potential file and initialize potential parameters
241
242 read_file(arg[2]);
243 setup_params();
244 }
245
246 /* ----------------------------------------------------------------------
247 init specific to this pair style
248 ------------------------------------------------------------------------- */
249
init_style()250 void PairAGNI::init_style()
251 {
252 // need a full neighbor list
253
254 int irequest = neighbor->request(this,instance_me);
255 neighbor->requests[irequest]->half = 0;
256 neighbor->requests[irequest]->full = 1;
257 }
258
259 /* ----------------------------------------------------------------------
260 init for one type pair i,j and corresponding j,i
261 ------------------------------------------------------------------------- */
262
init_one(int i,int j)263 double PairAGNI::init_one(int i, int j)
264 {
265 if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
266
267 return cutmax;
268 }
269
270 /* ---------------------------------------------------------------------- */
271
read_file(char * filename)272 void PairAGNI::read_file(char *filename)
273 {
274 int i,j,k,curparam,wantdata,fp_counter;
275
276 if (params) {
277 for (int i = 0; i < nparams; ++i) {
278 memory->destroy(params[i].eta);
279 memory->destroy(params[i].alpha);
280 memory->destroy(params[i].xU);
281 }
282 memory->destroy(params);
283 params = nullptr;
284 }
285 nparams = 0;
286
287 fp_counter = 0;
288 wantdata = -1;
289
290 // read potential file
291 if (comm->me == 0) {
292 PotentialFileReader reader(lmp, filename, "agni", unit_convert_flag);
293
294 char * line;
295
296 while ((line = reader.next_line())) {
297 try {
298 ValueTokenizer values(line);
299 if (wantdata == -1) {
300 std::string tag = values.next_string();
301
302 if (tag == "n_elements") {
303 nparams = values.next_int();
304
305 if ((nparams < 1) || params) // sanity check
306 error->all(FLERR,"Invalid AGNI potential file");
307 params = memory->create(params,nparams,"pair:params");
308 memset(params,0,nparams*sizeof(Param));
309
310 curparam = -1;
311 wantdata = -1;
312 } else if (tag == "interaction") {
313 for (j = 0; j < nparams; ++j) {
314 std::string element = values.next_string();
315 if (element == elements[params[j].ielement])
316 curparam = j;
317 }
318 } else if (tag == "element") {
319 for (j = 0; j < nparams; ++j) {
320 std::string element = values.next_string();
321 for (k = 0; k < nelements; ++k)
322 if (element == elements[k]) break;
323 if (k == nelements)
324 error->all(FLERR,"No suitable parameters for requested element found");
325 else params[j].ielement = k;
326 }
327 } else if (tag == "generation") {
328 atomic_feature_version = values.next_int();
329 if (atomic_feature_version != AGNI_VERSION_1 && atomic_feature_version != AGNI_VERSION_2)
330 error->all(FLERR,"Incompatible AGNI potential file version");
331 } else if (tag == "eta") {
332 params[curparam].numeta = values.count() - 1;
333 memory->create(params[curparam].eta,params[curparam].numeta,"agni:eta");
334 for (j = 0; j < params[curparam].numeta; j++)
335 params[curparam].eta[j] = values.next_double();
336 } else if (tag == "gwidth") {
337 params[curparam].gwidth = values.next_double();
338 } else if (tag == "Rc") {
339 params[curparam].cut = values.next_double();
340 } else if (tag == "n_train") {
341 params[curparam].numtrain = values.next_int();
342 memory->create(params[curparam].alpha,params[curparam].numtrain,"agni:alpha");
343 memory->create(params[curparam].xU,params[curparam].numtrain,params[curparam].numtrain,"agni:xU");
344 } else if (tag == "sigma") {
345 params[curparam].sigma = values.next_double();
346 } else if (tag == "b") {
347 params[curparam].b = values.next_double();
348 } else if (tag == "endVar") {
349 if (atomic_feature_version == AGNI_VERSION_1)
350 params[curparam].gwidth = 0.0;
351 wantdata = curparam;
352 curparam = -1;
353 } else error->warning(FLERR,"Ignoring unknown tag '{}' in AGNI potential file.",tag);
354 } else {
355 if (params && wantdata >=0) {
356 if ((int)values.count() == params[wantdata].numeta + 2) {
357 for (k = 0; k < params[wantdata].numeta; ++k)
358 params[wantdata].xU[k][fp_counter] = values.next_double();
359 values.next_double(); // ignore
360 params[wantdata].alpha[fp_counter] = values.next_double();
361 fp_counter++;
362 } else if ((int)values.count() == params[wantdata].numeta + 3) {
363 values.next_double(); // ignore
364 for (k = 0; k < params[wantdata].numeta; ++k)
365 params[wantdata].xU[k][fp_counter] = values.next_double();
366 values.next_double(); // ignore
367 params[wantdata].alpha[fp_counter] = values.next_double();
368 fp_counter++;
369 } else error->all(FLERR,"Invalid AGNI potential file");
370 }
371 }
372 }
373 catch (TokenizerException &e) {
374 error->one(FLERR, e.what());
375 }
376 }
377 }
378 MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
379 MPI_Bcast(&atomic_feature_version, 1, MPI_INT, 0, world);
380 if (comm->me != 0) {
381 params = memory->create(params,nparams,"pair:params");
382 memset(params,0,nparams*sizeof(Param));
383 }
384 MPI_Bcast(params, nparams*sizeof(Param), MPI_BYTE, 0, world);
385 for (i = 0; i < nparams; i++) {
386 if (comm->me != 0) {
387 memory->create(params[i].alpha,params[i].numtrain,"agni:alpha");
388 memory->create(params[i].eta,params[i].numeta,"agni:eta");
389 memory->create(params[i].xU,params[i].numeta,params[i].numtrain,"agni:xU");
390 }
391
392 MPI_Bcast(params[i].alpha, params[i].numtrain, MPI_DOUBLE, 0, world);
393 MPI_Bcast(params[i].eta, params[i].numeta, MPI_DOUBLE, 0, world);
394 MPI_Bcast(¶ms[i].xU[0][0],params[i].numtrain*params[i].numeta,MPI_DOUBLE,0,world);
395 }
396 }
397
398 /* ---------------------------------------------------------------------- */
399
setup_params()400 void PairAGNI::setup_params()
401 {
402 int i,m,n;
403 double rtmp;
404
405 // set elem1param for all elements
406
407 memory->destroy(elem1param);
408 memory->create(elem1param,nelements,"pair:elem1param");
409
410 for (i = 0; i < nelements; i++) {
411 n = -1;
412 for (m = 0; m < nparams; m++) {
413 if (i == params[m].ielement) {
414 if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
415 n = m;
416 }
417 }
418 if (n < 0) error->all(FLERR,"Potential file is missing an entry");
419 elem1param[i] = n;
420 }
421
422 // compute parameter values derived from inputs
423
424 // set cutsq using shortcut to reduce neighbor list for accelerated
425 // calculations. cut must remain unchanged as it is a potential parameter
426 // (cut = a*sigma)
427
428 cutmax = 0.0;
429 for (m = 0; m < nparams; m++) {
430 rtmp = params[m].cut;
431 params[m].cutsq = rtmp * rtmp;
432 if (rtmp > cutmax) cutmax = rtmp;
433 }
434 }
435