1 // clang-format off
2 /* -*- c++ -*- ----------------------------------------------------------
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 Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
16 Doyl Dickel (MSU) doyl@me.msstate.edu
17 ----------------------------------------------------------------------*/
18 /*
19 “The research described and the resulting data presented herein, unless
20 otherwise noted, was funded under PE 0602784A, Project T53 "Military
21 Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
22 managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
23 the Engineer Research and Development Center (ERDC). The work described in
24 this document was conducted at CAVS, MSU. Permission was granted by ERDC
25 to publish this information. Any opinions, findings and conclusions or
26 recommendations expressed in this material are those of the author(s) and
27 do not necessarily reflect the views of the United States Army.”
28
29 DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
30 */
31
32 #include "pair_rann.h"
33
34 #include "atom.h"
35 #include "citeme.h"
36 #include "error.h"
37 #include "math_special.h"
38 #include "memory.h"
39 #include "neighbor.h"
40 #include "neigh_list.h"
41 #include "neigh_request.h"
42 #include "tokenizer.h"
43 #include "update.h"
44
45 #include <cmath>
46 #include <cstring>
47
48 #include "rann_activation_linear.h"
49 #include "rann_activation_sig_i.h"
50 #include "rann_fingerprint_bond.h"
51 #include "rann_fingerprint_bondscreened.h"
52 #include "rann_fingerprint_bondscreenedspin.h"
53 #include "rann_fingerprint_bondspin.h"
54 #include "rann_fingerprint_radial.h"
55 #include "rann_fingerprint_radialscreened.h"
56 #include "rann_fingerprint_radialscreenedspin.h"
57 #include "rann_fingerprint_radialspin.h"
58
59 #define MAXLINE 1024
60
61 using namespace LAMMPS_NS;
62
63 static const char cite_ml_rann_package[] =
64 "ML-RANN package:\n\n"
65 "@Article{Nitol2021,\n"
66 " author = {Nitol, Mashroor S and Dickel, Doyl E and Barrett, Christopher D},\n"
67 " title = {Artificial neural network potential for pure zinc},\n"
68 " journal = {Computational Materials Science},\n"
69 " year = 2021,\n"
70 " volume = 188,\n"
71 " pages = {110207}\n"
72 "}\n\n";
73
74
PairRANN(LAMMPS * lmp)75 PairRANN::PairRANN(LAMMPS *lmp) : Pair(lmp)
76 {
77 if (lmp->citeme) lmp->citeme->add(cite_ml_rann_package);
78
79 //initialize ints and bools
80 single_enable = 0;
81 restartinfo = 0;
82 one_coeff = 1;
83 manybody_flag = 1;
84 allocated = 0;
85 nelements = -1;
86 nelementsp = -1;
87 comm_forward = 0;
88 comm_reverse = 0;
89 res = 10000;
90 cutmax = 0;
91 dospin = false;
92 memguess = 0;
93 nmax1 = 0;
94 nmax2 = 0;
95 fmax = 0;
96 fnmax = 0;
97 //at least one of the following two flags will change during fingerprint definition:
98 doscreen = false;
99 allscreen = true;
100
101 //null init for arrays with sizes not yet determined.
102 elements = nullptr;
103 mass = nullptr;
104 elementsp = nullptr;
105 map = nullptr;
106 fingerprintcount = nullptr;
107 fingerprintlength = nullptr;
108 fingerprintperelement = nullptr;
109 screening_min = nullptr;
110 screening_max = nullptr;
111 weightdefined = nullptr;
112 biasdefined = nullptr;
113 xn = nullptr;
114 yn = nullptr;
115 zn = nullptr;
116 Sik = nullptr;
117 dSikx = nullptr;
118 dSiky = nullptr;
119 dSikz = nullptr;
120 dSijkx = nullptr;
121 dSijky = nullptr;
122 dSijkz = nullptr;
123 sx = nullptr;
124 sy = nullptr;
125 sz = nullptr;
126 dSijkxc = nullptr;
127 dSijkyc = nullptr;
128 dSijkzc = nullptr;
129 dfeaturesx = nullptr;
130 dfeaturesy = nullptr;
131 dfeaturesz = nullptr;
132 features = nullptr;
133 layer = nullptr;
134 sum = nullptr;
135 sum1 = nullptr;
136 dlayerx = nullptr;
137 dlayery = nullptr;
138 dlayerz = nullptr;
139 dlayersumx = nullptr;
140 dlayersumy = nullptr;
141 dlayersumz = nullptr;
142 dsx = nullptr;
143 dsy = nullptr;
144 dsz = nullptr;
145 dssumx = nullptr;
146 dssumy = nullptr;
147 dssumz = nullptr;
148 tn = nullptr;
149 jl = nullptr;
150 Bij = nullptr;
151 sims = nullptr;
152 net = nullptr;
153 activation = nullptr;
154 fingerprints = nullptr;
155 }
156
~PairRANN()157 PairRANN::~PairRANN()
158 {
159 deallocate();
160 }
161
deallocate()162 void PairRANN::deallocate()
163 {
164 //clear memory
165 delete[] mass;
166 for (int i=0;i<nelements;i++) {delete [] elements[i];}
167 delete[] elements;
168 for (int i=0;i<nelementsp;i++) {delete [] elementsp[i];}
169 delete[] elementsp;
170 for (int i=0;i<=nelements;i++) {
171 if (net[i].layers>0) {
172 for (int j=0;j<net[i].layers-1;j++) {
173 delete[] net[i].Weights[j];
174 delete[] net[i].Biases[j];
175 delete activation[i][j];
176 }
177 delete[] activation[i];
178 delete[] net[i].Weights;
179 delete[] net[i].Biases;
180 delete[] weightdefined[i];
181 delete[] biasdefined[i];
182 }
183 delete[] net[i].dimensions;
184 }
185 delete[] net;
186 delete[] map;
187 for (int i=0;i<nelementsp;i++) {
188 if (fingerprintlength[i]>0) {
189 for (int j=0;j<fingerprintperelement[i];j++) {
190 delete fingerprints[i][j];
191 }
192 delete[] fingerprints[i];
193 }
194 }
195 delete[] fingerprints;
196 delete[] activation;
197 delete[] weightdefined;
198 delete[] biasdefined;
199 delete[] fingerprintcount;
200 delete[] fingerprintperelement;
201 delete[] fingerprintlength;
202 delete[] screening_min;
203 delete[] screening_max;
204 memory->destroy(xn);
205 memory->destroy(yn);
206 memory->destroy(zn);
207 memory->destroy(tn);
208 memory->destroy(jl);
209 memory->destroy(features);
210 memory->destroy(dfeaturesx);
211 memory->destroy(dfeaturesy);
212 memory->destroy(dfeaturesz);
213 memory->destroy(layer);
214 memory->destroy(sum);
215 memory->destroy(sum1);
216 memory->destroy(dlayerx);
217 memory->destroy(dlayery);
218 memory->destroy(dlayerz);
219 memory->destroy(dlayersumx);
220 memory->destroy(dlayersumy);
221 memory->destroy(dlayersumz);
222 memory->destroy(Sik);
223 memory->destroy(Bij);
224 memory->destroy(dSikx);
225 memory->destroy(dSiky);
226 memory->destroy(dSikz);
227 memory->destroy(dSijkx);
228 memory->destroy(dSijky);
229 memory->destroy(dSijkz);
230 memory->destroy(dSijkxc);
231 memory->destroy(dSijkyc);
232 memory->destroy(dSijkzc);
233 memory->destroy(sx);
234 memory->destroy(sy);
235 memory->destroy(sz);
236 memory->destroy(dsx);
237 memory->destroy(dsy);
238 memory->destroy(dsz);
239 memory->destroy(dssumx);
240 memory->destroy(dssumy);
241 memory->destroy(dssumz);
242 memory->destroy(setflag);
243 memory->destroy(cutsq);
244 }
245
allocate(const std::vector<std::string> & elementwords)246 void PairRANN::allocate(const std::vector<std::string> &elementwords)
247 {
248 int i,n;
249 n = atom->ntypes;
250 memory->create(setflag,n+1,n+1,"pair:setflag");
251 memory->create(cutsq,n+1,n+1,"pair:cutsq");
252 cutmax = 0;
253 nmax1 = 100;
254 nmax2 = 20;
255 fmax = 0;
256 fnmax = 0;
257 nelementsp=nelements+1;
258 //initialize arrays
259 elements = new char *[nelements];
260 elementsp = new char *[nelementsp];//elements + 'all'
261 mass = new double[nelements];
262 net = new NNarchitecture[nelementsp];
263 weightdefined = new bool*[nelementsp];
264 biasdefined = new bool *[nelementsp];
265 activation = new RANN::Activation**[nelementsp];
266 fingerprints = new RANN::Fingerprint**[nelementsp];
267 fingerprintlength = new int[nelementsp];
268 fingerprintperelement = new int[nelementsp];
269 fingerprintcount = new int[nelementsp];
270 screening_min = new double[nelements*nelements*nelements];
271 screening_max = new double[nelements*nelements*nelements];
272 for (i=0;i<nelements;i++) {
273 for (int j =0;j<nelements;j++) {
274 for (int k=0;k<nelements;k++) {
275 screening_min[i*nelements*nelements+j*nelements+k] = 0.8;//default values. Custom values may be read from potential file later.
276 screening_max[i*nelements*nelements+j*nelements+k] = 2.8;//default values. Custom values may be read from potential file later.
277 }
278 }
279 }
280 for (i=0;i<=nelements;i++) {
281 fingerprintlength[i]=0;
282 fingerprintperelement[i] = -1;
283 fingerprintcount[i] = 0;
284 if (i<nelements) {
285 mass[i]=-1.0;
286 elements[i]= utils::strdup(elementwords[i]);
287 }
288 elementsp[i] = utils::strdup(elementwords[i]);
289 net[i].layers = 0;
290 net[i].dimensions = new int[1];
291 net[i].dimensions[0]=0;
292 }
293 }
294
settings(int narg,char **)295 void PairRANN::settings(int narg, char ** /*arg*/)
296 {
297 //read pair_style command in input file
298 if (narg > 0) error->one(FLERR,"Illegal pair_style command");
299 }
300
coeff(int narg,char ** arg)301 void PairRANN::coeff(int narg, char **arg)
302 {
303 int i,j;
304 deallocate();//clear allocation from any previous coeff
305 map = new int[atom->ntypes+1];
306 if (narg != 3 + atom->ntypes) error->one(FLERR,"Incorrect args for pair coefficients");
307 if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) error->one(FLERR,"Incorrect args for pair coefficients");
308 nelements = -1;
309 read_file(arg[2]);
310 // read args that map atom types to elements in potential file
311 // map[i] = which element the Ith atom type is, -1 if NULL
312 for (i = 3; i < narg; i++) {
313 if (strcmp(arg[i],"NULL") == 0) {
314 map[i-2] = -1;
315 continue;
316 }
317 for (j = 0; j < nelements; j++) {
318 if (strcmp(arg[i],elements[j]) == 0) break;
319 }
320 if (j < nelements) map[i-2] = j;
321 else error->one(FLERR,"No matching element in NN potential file");
322 }
323 // clear setflag since coeff() called once with I,J = * *
324 int n = atom->ntypes;
325 for (i = 1; i <= n; i++) {
326 for (j = i; j <= n; j++) {
327 setflag[i][j] = 0;
328 }
329 }
330 // set setflag i,j for type pairs where both are mapped to elements
331 // set mass of atom type if i = j
332 int count = 0;
333 for (i = 1; i <= n; i++) {
334 for (j = i; j <= n; j++) {
335 if (map[i] >= 0 && map[j] >= 0) {
336 setflag[i][j] = 1;
337 if (i == j) atom->set_mass(FLERR,i,mass[map[i]]);
338 count++;
339 }
340 }
341 }
342 if (count == 0) error->one(FLERR,"Incorrect args for pair coefficients");
343 for (i=0;i<nelementsp;i++) {
344 for (j=0;j<fingerprintperelement[i];j++) {
345 fingerprints[i][j]->allocate();
346 }
347 }
348 allocated=1;
349 }
350
read_file(char * filename)351 void PairRANN::read_file(char *filename)
352 {
353 FILE *fp;
354 int eof = 0;
355 std::string line,line1;
356 const int longline = 4096;
357 int linenum=0;
358 char linetemp[longline];
359 std::string strtemp;
360 char *ptr;
361 std::vector<std::string> linev,line1v;
362 fp = utils::open_potential(filename,lmp,nullptr);
363 if (fp == nullptr) {error->one(FLERR,"Cannot open RANN potential file");}
364 ptr=fgets(linetemp,longline,fp);
365 linenum++;
366 strtemp=utils::trim_comment(linetemp);
367 while (strtemp.empty()) {
368 ptr=fgets(linetemp,longline,fp);
369 strtemp=utils::trim_comment(linetemp);
370 linenum++;
371 }
372 line=strtemp;
373 while (eof == 0) {
374 ptr=fgets(linetemp,longline,fp);
375 linenum++;
376 if (ptr == nullptr) {
377 fclose(fp);
378 if (check_potential()) {
379 error->one(FLERR,"Invalid syntax in potential file, values are inconsistent or missing");
380 }
381 else {
382 eof = 1;
383 break;
384 }
385 }
386 strtemp=utils::trim_comment(linetemp);
387 while (strtemp.empty()) {
388 ptr=fgets(linetemp,longline,fp);
389 strtemp=utils::trim_comment(linetemp);
390 linenum++;
391 }
392 line1=linetemp;
393 Tokenizer values = Tokenizer(line,": ,\t_\n");
394 Tokenizer values1 = Tokenizer(line1,": ,\t_\n");
395 linev = values.as_vector();
396 line1v = values1.as_vector();
397 if (linev[0]=="atomtypes") read_atom_types(line1v,filename,linenum);
398 else if (linev[0]=="mass") read_mass(linev,line1v,filename,linenum);
399 else if (linev[0]=="fingerprintsperelement") read_fpe(linev,line1v,filename,linenum);
400 else if (linev[0]=="fingerprints") read_fingerprints(linev,line1v,filename,linenum);
401 else if (linev[0]=="fingerprintconstants") read_fingerprint_constants(linev,line1v,filename,linenum);
402 else if (linev[0]=="networklayers") read_network_layers(linev,line1v,filename,linenum);
403 else if (linev[0]=="layersize") read_layer_size(linev,line1v,filename,linenum);
404 else if (linev[0]=="weight") read_weight(linev,line1v,fp,filename,&linenum);
405 else if (linev[0]=="bias") read_bias(linev,line1v,fp,filename,&linenum);
406 else if (linev[0]=="activationfunctions") read_activation_functions(linev,line1v,filename,linenum);
407 else if (linev[0]=="screening") read_screening(linev,line1v,filename,linenum);
408 else if (linev[0]=="calibrationparameters") continue;//information on how the network was trained
409 else error->one(FLERR,"Could not understand file syntax: unknown keyword");
410 ptr=fgets(linetemp,longline,fp);
411 linenum++;
412 strtemp=utils::trim_comment(linetemp);
413 while (strtemp.empty()) {
414 ptr=fgets(linetemp,longline,fp);
415 strtemp=utils::trim_comment(linetemp);
416 linenum++;
417 }
418 if (ptr == nullptr) {
419 if (check_potential()) {
420 error->one(FLERR,"Invalid syntax in potential file, values are inconsistent or missing");
421 }
422 else {
423 eof = 1;
424 break;
425 }
426 }
427 line=linetemp;
428 }
429 }
430
read_atom_types(std::vector<std::string> line,char * filename,int linenum)431 void PairRANN::read_atom_types(std::vector<std::string> line,char *filename,int linenum) {
432 int nwords = line.size();
433 if (nwords < 1) error->one(filename,linenum,"Incorrect syntax for atom types");
434 nelements = nwords;
435 line.push_back("all");
436 allocate(line);
437 }
438
read_mass(const std::vector<std::string> & line1,const std::vector<std::string> & line2,const char * filename,int linenum)439 void PairRANN::read_mass(const std::vector<std::string> &line1, const std::vector<std::string> &line2, const char *filename,int linenum) {
440 if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before mass in potential file.");
441 for (int i=0;i<nelements;i++) {
442 if (line1[1].compare(elements[i])==0) {
443 mass[i]=utils::numeric(filename,linenum,line2[0].c_str(),1,lmp);
444 return;
445 }
446 }
447 error->one(filename,linenum-1,"mass element not found in atom types.");
448 }
449
read_fpe(std::vector<std::string> line,std::vector<std::string> line1,char * filename,int linenum)450 void PairRANN::read_fpe(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
451 int i;
452 if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before fingerprints per element in potential file.");
453 for (i=0;i<nelementsp;i++) {
454 if (line[1].compare(elementsp[i])==0) {
455 fingerprintperelement[i] = utils::inumeric(filename,linenum,line1[0].c_str(),1,lmp);
456 fingerprints[i] = new RANN::Fingerprint *[fingerprintperelement[i]];
457 for (int j=0;j<fingerprintperelement[i];j++) {
458 fingerprints[i][j]=new RANN::Fingerprint(this);
459 }
460 return;
461 }
462 }
463 error->one(filename,linenum-1,"fingerprint-per-element element not found in atom types");
464 }
465
read_fingerprints(std::vector<std::string> line,std::vector<std::string> line1,char * filename,int linenum)466 void PairRANN::read_fingerprints(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
467 int nwords1,nwords,i,j,k,i1,*atomtypes;
468 bool found;
469 nwords1 = line1.size();
470 nwords = line.size();
471 if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before fingerprints in potential file.");
472 atomtypes = new int[nwords-1];
473 for (i=1;i<nwords;i++) {
474 found = false;
475 for (j=0;j<nelementsp;j++) {
476 if (line[i].compare(elementsp[j])==0) {
477 atomtypes[i-1]=j;
478 found = true;
479 break;
480 }
481 }
482 if (!found) {error->one(filename,linenum-1,"fingerprint element not found in atom types");}
483 }
484 i = atomtypes[0];
485 k = 0;
486 if (fingerprintperelement[i]==-1) {error->one(filename,linenum-1,"fingerprint per element must be defined before fingerprints");}
487 while (k<nwords1) {
488 i1 = fingerprintcount[i];
489 if (i1>=fingerprintperelement[i]) {error->one(filename,linenum,"more fingerprints defined than fingerprint per element");}
490 delete fingerprints[i][i1];
491 fingerprints[i][i1] = create_fingerprint(line1[k].c_str());
492 if (fingerprints[i][i1]->n_body_type!=nwords-1) {error->one(filename,linenum,"invalid fingerprint for element combination");}
493 k++;
494 fingerprints[i][i1]->init(atomtypes,utils::inumeric(filename,linenum,line1[k++].c_str(),1,lmp));
495 fingerprintcount[i]++;
496 }
497 delete[] atomtypes;
498 }
499
read_fingerprint_constants(std::vector<std::string> line,std::vector<std::string> line1,char * filename,int linenum)500 void PairRANN::read_fingerprint_constants(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
501 int i,j,k,i1,*atomtypes;
502 bool found;
503 int nwords = line.size();
504 if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before fingerprints in potential file.");
505 int n_body_type = nwords-4;
506 atomtypes = new int[n_body_type];
507 for (i=1;i<=n_body_type;i++) {
508 found = false;
509 for (j=0;j<nelementsp;j++) {
510 if (line[i].compare(elementsp[j])==0) {
511 atomtypes[i-1]=j;
512 found = true;
513 break;
514 }
515 }
516 if (!found) {error->one(filename,linenum-1,"fingerprint element not found in atom types");}
517 }
518 i = atomtypes[0];
519 found = false;
520 for (k=0;k<fingerprintperelement[i];k++) {
521 if (fingerprints[i][k]->empty) {continue;}
522 if (n_body_type!=fingerprints[i][k]->n_body_type) {continue;}
523 for (j=0;j<n_body_type;j++) {
524 if (fingerprints[i][k]->atomtypes[j]!=atomtypes[j]) {break;}
525 if (j==n_body_type-1) {
526 if (line[nwords-3].compare(fingerprints[i][k]->style)==0 && utils::inumeric(filename,linenum,line[nwords-2].c_str(),1,lmp)==fingerprints[i][k]->id) {
527 found=true;
528 i1 = k;
529 break;
530 }
531 }
532 }
533 if (found) {break;}
534 }
535 if (!found) {error->one(filename,linenum-1,"cannot define constants for unknown fingerprint");}
536 fingerprints[i][i1]->fullydefined=fingerprints[i][i1]->parse_values(line[nwords-1],line1);
537 delete[] atomtypes;
538 }
539
read_network_layers(std::vector<std::string> line,std::vector<std::string> line1,char * filename,int linenum)540 void PairRANN::read_network_layers(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
541 int i,j;
542 if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before network layers in potential file.");
543 for (i=0;i<nelements;i++) {
544 if (line[1].compare(elements[i])==0) {
545 net[i].layers = utils::inumeric(filename,linenum,line1[0].c_str(),1,lmp);
546 if (net[i].layers < 1)error->one(filename,linenum,"invalid number of network layers");
547 delete[] net[i].dimensions;
548 weightdefined[i] = new bool [net[i].layers];
549 biasdefined[i] = new bool [net[i].layers];
550 net[i].dimensions = new int[net[i].layers];
551 net[i].Weights = new double * [net[i].layers-1];
552 net[i].Biases = new double * [net[i].layers-1];
553 for (j=0;j<net[i].layers;j++) {
554 net[i].dimensions[j]=0;
555 weightdefined[i][j] = false;
556 biasdefined[i][j] = false;
557 }
558 activation[i]=new RANN::Activation* [net[i].layers-1];
559 for (int j=0;j<net[i].layers-1;j++) {
560 activation[i][j]= new RANN::Activation(this);
561 }
562 return;
563 }
564 }
565 error->one(filename,linenum-1,"network layers element not found in atom types");
566 }
567
read_layer_size(std::vector<std::string> line,std::vector<std::string> line1,char * filename,int linenum)568 void PairRANN::read_layer_size(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
569 int i;
570 for (i=0;i<nelements;i++) {
571 if (line[1].compare(elements[i])==0) {
572 if (net[i].layers==0)error->one(filename,linenum-1,"networklayers for each atom type must be defined before the corresponding layer sizes.");
573 int j = utils::inumeric(filename,linenum,line[2].c_str(),1,lmp);
574 if (j>=net[i].layers || j<0) {error->one(filename,linenum,"invalid layer in layer size definition");};
575 net[i].dimensions[j]= utils::inumeric(filename,linenum,line1[0].c_str(),1,lmp);
576 return;
577 }
578 }
579 error->one(filename,linenum-1,"layer size element not found in atom types");
580 }
581
read_weight(std::vector<std::string> line,std::vector<std::string> line1,FILE * fp,char * filename,int * linenum)582 void PairRANN::read_weight(std::vector<std::string> line,std::vector<std::string> line1,FILE* fp,char *filename,int *linenum) {
583 int i,j,k,l,nwords;
584 char *ptr;
585 const int longline = 4096;
586 char linetemp [longline];
587 for (l=0;l<nelements;l++) {
588 if (line[1].compare(elements[l])==0) {
589 if (net[l].layers==0)error->one(filename,*linenum-1,"networklayers must be defined before weights.");
590 i=utils::inumeric(filename,*linenum,line[2].c_str(),1,lmp);
591 if (i>=net[l].layers || i<0)error->one(filename,*linenum-1,"invalid weight layer");
592 if (net[l].dimensions[i]==0 || net[l].dimensions[i+1]==0) error->one(filename,*linenum-1,"network layer sizes must be defined before corresponding weight");
593 net[l].Weights[i] = new double[net[l].dimensions[i]*net[l].dimensions[i+1]];
594 weightdefined[l][i] = true;
595 nwords = line1.size();
596 if (nwords != net[l].dimensions[i])error->one(filename,*linenum,"invalid weights per line");
597 for (k=0;k<net[l].dimensions[i];k++) {
598 net[l].Weights[i][k] = utils::numeric(filename,*linenum,line1[k].c_str(),1,lmp);
599 }
600 for (j=1;j<net[l].dimensions[i+1];j++) {
601 ptr = fgets(linetemp,longline,fp);
602 (*linenum)++;
603 Tokenizer values1 = Tokenizer(linetemp,": ,\t_\n");
604 line1 = values1.as_vector();
605 if (ptr==nullptr)error->one(filename,*linenum,"unexpected end of potential file!");
606 nwords = line1.size();
607 if (nwords != net[l].dimensions[i])error->one(filename,*linenum,"invalid weights per line");
608 for (k=0;k<net[l].dimensions[i];k++) {
609 net[l].Weights[i][j*net[l].dimensions[i]+k] = utils::numeric(filename,*linenum,line1[k].c_str(),1,lmp);
610 }
611 }
612 return;
613 }
614 }
615 error->one(filename,*linenum-1,"weight element not found in atom types");
616 }
617
read_bias(std::vector<std::string> line,std::vector<std::string> line1,FILE * fp,char * filename,int * linenum)618 void PairRANN::read_bias(std::vector<std::string> line,std::vector<std::string> line1,FILE* fp,char *filename,int *linenum) {
619 int i,j,l;
620 char linetemp[MAXLINE],*ptr;
621 for (l=0;l<nelements;l++) {
622 if (line[1].compare(elements[l])==0) {
623 if (net[l].layers==0)error->one(filename,*linenum-1,"networklayers must be defined before biases.");
624 i=utils::inumeric(filename,*linenum,line[2].c_str(),1,lmp);
625 if (i>=net[l].layers || i<0)error->one(filename,*linenum-1,"invalid bias layer");
626 if (net[l].dimensions[i]==0) error->one(filename,*linenum-1,"network layer sizes must be defined before corresponding bias");
627 biasdefined[l][i] = true;
628 net[l].Biases[i] = new double[net[l].dimensions[i+1]];
629 net[l].Biases[i][0] = utils::numeric(filename,*linenum,line1[0].c_str(),1,lmp);
630 for (j=1;j<net[l].dimensions[i+1];j++) {
631 ptr=fgets(linetemp,MAXLINE,fp);
632 if (ptr==nullptr)error->one(filename,*linenum,"unexpected end of potential file!");
633 (*linenum)++;
634 Tokenizer values1 = Tokenizer(linetemp,": ,\t_\n");
635 line1 = values1.as_vector();
636 net[l].Biases[i][j] = utils::numeric(filename,*linenum,line1[0].c_str(),1,lmp);
637 }
638 return;
639 }
640 }
641 error->one(filename,*linenum-1,"bias element not found in atom types");
642 }
643
read_activation_functions(std::vector<std::string> line,std::vector<std::string> line1,char * filename,int linenum)644 void PairRANN::read_activation_functions(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
645 int i,l;
646 for (l=0;l<nelements;l++) {
647 if (line[1].compare(elements[l])==0) {
648 if (net[l].layers==0)error->one(filename,linenum-1,"networklayers must be defined before activation functions.");
649 i = strtol(line[2].c_str(),nullptr,10);
650 if (i>=net[l].layers || i<0)error->one(filename,linenum-1,"invalid activation layer");
651 delete activation[l][i];
652 activation[l][i]=create_activation(line1[0].c_str());
653 return;
654 }
655 }
656 error->one(filename,linenum-1,"activation function element not found in atom types");
657 }
658
read_screening(std::vector<std::string> line,std::vector<std::string> line1,char * filename,int linenum)659 void PairRANN::read_screening(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
660 int i,j,k,*atomtypes;
661 bool found;
662 int nwords = line.size();
663 if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before fingerprints in potential file.");
664 if (nwords!=5)error->one(filename,linenum-1,"invalid screening command");
665 int n_body_type = 3;
666 atomtypes = new int[n_body_type];
667 for (i=1;i<=n_body_type;i++) {
668 found = false;
669 for (j=0;j<nelementsp;j++) {
670 if (line[i].compare(elementsp[j])==0) {
671 atomtypes[i-1]=j;
672 found = true;
673 break;
674 }
675 }
676 if (!found) {error->one(filename,linenum-1,"fingerprint element not found in atom types");}
677 }
678 i = atomtypes[0];
679 j = atomtypes[1];
680 k = atomtypes[2];
681 int index = i*nelements*nelements+j*nelements+k;
682 if (line[4].compare("Cmin")==0) {
683 screening_min[index] = utils::numeric(filename,linenum,line1[0].c_str(),1,lmp);
684 }
685 else if (line[4].compare("Cmax")==0) {
686 screening_max[index] = utils::numeric(filename,linenum,line1[0].c_str(),1,lmp);
687 }
688 else error->one(filename,linenum-1,"unrecognized screening keyword");
689 delete[] atomtypes;
690 }
691
692 //Called after finishing reading the potential file to make sure it is complete. True is bad.
693 //also does the rest of the memory allocation.
check_potential()694 bool PairRANN::check_potential() {
695 int i,j,k,l;
696 if (nelements==-1) {return true;}
697 for (i=0;i<=nelements;i++) {
698 if (i<nelements) {
699 if (mass[i]<0)return true;//uninitialized mass
700 }
701 if (net[i].layers==0)break;//no definitions for this starting element, not considered an error.
702 net[i].maxlayer=0;
703 for (j=0;j<net[i].layers;j++) {
704 if (net[i].dimensions[j]==0)return true;//incomplete network definition
705 if (net[i].dimensions[j]>net[i].maxlayer)net[i].maxlayer = net[i].dimensions[j];
706 }
707 if (net[i].maxlayer>fnmax) {fnmax = net[i].maxlayer;}
708 if (net[i].dimensions[net[i].layers-1]!=1)return true;//output layer must have single neuron (the energy)
709 if (net[i].dimensions[0]>fmax)fmax=net[i].dimensions[0];
710 for (j=0;j<net[i].layers-1;j++) {
711 if (!weightdefined[i][j])return true;//undefined weights
712 if (!biasdefined[i][j])return true;//undefined biases
713 if (activation[i][j]->empty)return true;//undefined activations
714 for (k=0;k<net[i].dimensions[j+1];k++) {
715 for (l=0;l<net[i].dimensions[j];l++) {
716 if (net[i].Weights[j][k*net[i].dimensions[j]+l]==0)return true;//undefined weights
717 }
718 if (net[i].Biases[j][k]==0)return true;//undefined biases
719 }
720 }
721 for (j=0;j<fingerprintperelement[i];j++) {
722 if (fingerprints[i][j]->fullydefined==false)return true;
723 fingerprints[i][j]->startingneuron = fingerprintlength[i];
724 fingerprintlength[i] +=fingerprints[i][j]->get_length();
725 if (fingerprints[i][j]->rc>cutmax) {cutmax = fingerprints[i][j]->rc;}
726 }
727 if (net[i].dimensions[0]!=fingerprintlength[i])return true;
728 }
729 memory->create(xn,nmax1,"pair:xn");
730 memory->create(yn,nmax1,"pair:yn");
731 memory->create(zn,nmax1,"pair:zn");
732 memory->create(tn,nmax1,"pair:tn");
733 memory->create(jl,nmax1,"pair:jl");
734 memory->create(features,fmax,"pair:features");
735 memory->create(dfeaturesx,fmax*nmax2,"pair:dfeaturesx");
736 memory->create(dfeaturesy,fmax*nmax2,"pair:dfeaturesy");
737 memory->create(dfeaturesz,fmax*nmax2,"pair:dfeaturesz");
738 memory->create(layer,fnmax,"pair:layer");
739 memory->create(sum,fnmax,"pair:sum");
740 memory->create(sum1,fnmax,"pair:sum1");
741 memory->create(dlayerx,nmax2,fnmax,"pair:dlayerx");
742 memory->create(dlayery,nmax2,fnmax,"pair:dlayery");
743 memory->create(dlayerz,nmax2,fnmax,"pair:dlayerz");
744 memory->create(dlayersumx,nmax2,fnmax,"pair:dlayersumx");
745 memory->create(dlayersumy,nmax2,fnmax,"pair:dlayersumy");
746 memory->create(dlayersumz,nmax2,fnmax,"pair:dlayersumz");
747 if (doscreen) {
748 memory->create(Sik,nmax2,"pair:Sik");
749 memory->create(Bij,nmax2,"pair:Bij");
750 memory->create(dSikx,nmax2,"pair:dSikx");
751 memory->create(dSiky,nmax2,"pair:dSiky");
752 memory->create(dSikz,nmax2,"pair:dSikz");
753 memory->create(dSijkx,nmax2*nmax2,"pair:dSijkx");
754 memory->create(dSijky,nmax2*nmax2,"pair:dSijky");
755 memory->create(dSijkz,nmax2*nmax2,"pair:dSijkz");
756 memory->create(dSijkxc,nmax2,nmax2,"pair:dSijkxc");
757 memory->create(dSijkyc,nmax2,nmax2,"pair:dSijkyc");
758 memory->create(dSijkzc,nmax2,nmax2,"pair:dSijkzc");
759 }
760 if (dospin) {
761 memory->create(sx,fmax*nmax2,"pair:sx");
762 memory->create(sy,fmax*nmax2,"pair:sy");
763 memory->create(sz,fmax*nmax2,"pair:sz");
764 memory->create(dsx,nmax2,fnmax,"pair:dsx");
765 memory->create(dsy,nmax2,fnmax,"pair:dsy");
766 memory->create(dsz,nmax2,fnmax,"pair:dsz");
767 memory->create(dssumx,nmax2,fnmax,"pair:dssumx");
768 memory->create(dssumy,nmax2,fnmax,"pair:dssumy");
769 memory->create(dssumz,nmax2,fnmax,"pair:dssumz");
770 }
771 return false;//everything looks good
772 }
773
compute(int eflag,int vflag)774 void PairRANN::compute(int eflag, int vflag)
775 {
776 //perform force/energy computation_
777 if (dospin) {
778 if (strcmp(update->unit_style,"metal") != 0)
779 error->one(FLERR,"Spin pair styles require metal units");
780 if (!atom->sp_flag)
781 error->one(FLERR,"Spin pair styles requires atom/spin style");
782 }
783 ev_init(eflag,vflag);
784
785 // only global virial via fdotr is supported by this pair style
786
787 if (vflag_atom)
788 error->all(FLERR,"Pair style rann does not support computing per-atom stress");
789 if (vflag && !vflag_fdotr)
790 error->all(FLERR,"Pair style rann does not support 'pair_modify nofdotr'");
791
792 int ii,i,j;
793 int nn = 0;
794 sims = new Simulation[1];
795 sims->inum = listfull->inum;
796 sims->ilist=listfull->ilist;
797 sims->id = listfull->ilist;
798 sims->type = atom->type;
799 sims->x = atom->x;
800 sims->numneigh = listfull->numneigh;
801 sims->firstneigh = listfull->firstneigh;
802 if (dospin) {
803 sims->s = atom->sp;
804 }
805 int itype,f,jnum,len;
806 if (eflag_global) {eng_vdwl=0;eng_coul=0;}
807 double energy=0;
808 double **force = atom->f;
809 double **fm = atom->fm;
810
811 //loop over atoms
812 for (ii=0;ii<sims->inum;ii++) {
813 i = sims->ilist[ii];
814 itype = map[sims->type[i]];
815 f = net[itype].dimensions[0];
816 jnum = sims->numneigh[i];
817 if (jnum>nmax1) {
818 nmax1 = jnum;
819 memory->grow(xn,nmax1,"pair:xn");
820 memory->grow(yn,nmax1,"pair:yn");
821 memory->grow(zn,nmax1,"pair:zn");
822 memory->grow(tn,nmax1,"pair:tn");
823 memory->grow(jl,nmax1,"pair:jl");
824 }
825 cull_neighbor_list(&jnum,i,0);
826 if (jnum>nmax2) {
827 nmax2=jnum;
828 memory->grow(dfeaturesx,fmax*nmax2,"pair:dfeaturesx");
829 memory->grow(dfeaturesy,fmax*nmax2,"pair:dfeaturesy");
830 memory->grow(dfeaturesz,fmax*nmax2,"pair:dfeaturesz");
831 memory->grow(layer,fnmax,"pair:layer");
832 memory->grow(sum,fnmax,"pair:sum");
833 memory->grow(sum1,fnmax,"pair:sum1");
834 memory->grow(dlayerx,nmax2,fnmax,"pair:dlayerx");
835 memory->grow(dlayery,nmax2,fnmax,"pair:dlayery");
836 memory->grow(dlayerz,nmax2,fnmax,"pair:dlayerz");
837 memory->grow(dlayersumx,nmax2,fnmax,"pair:dlayersumx");
838 memory->grow(dlayersumy,nmax2,fnmax,"pair:dlayersumy");
839 memory->grow(dlayersumz,nmax2,fnmax,"pair:dlayersumz");
840 if (doscreen) {
841 memory->grow(Sik,nmax2,"pair:Sik");
842 memory->grow(Bij,nmax2,"pair:Bij");
843 memory->grow(dSikx,nmax2,"pair:dSikx");
844 memory->grow(dSiky,nmax2,"pair:dSiky");
845 memory->grow(dSikz,nmax2,"pair:dSikz");
846 memory->grow(dSijkx,nmax2*nmax2,"pair:dSijkx");
847 memory->grow(dSijky,nmax2*nmax2,"pair:dSijky");
848 memory->grow(dSijkz,nmax2*nmax2,"pair:dSijkz");
849 memory->destroy(dSijkxc);
850 memory->destroy(dSijkyc);
851 memory->destroy(dSijkzc);
852 memory->create(dSijkxc,nmax2,nmax2,"pair:dSijkxc");
853 memory->create(dSijkyc,nmax2,nmax2,"pair:dSijkyc");
854 memory->create(dSijkzc,nmax2,nmax2,"pair:dSijkzc");
855 }
856 if (dospin) {
857 memory->grow(sx,fmax*nmax2,"pair:sx");
858 memory->grow(sy,fmax*nmax2,"pair:sy");
859 memory->grow(sz,fmax*nmax2,"pair:sz");
860 memory->grow(dsx,nmax2,fnmax,"pair:dsx");
861 memory->grow(dsy,nmax2,fnmax,"pair:dsy");
862 memory->grow(dsz,nmax2,fnmax,"pair:dsz");
863 memory->grow(dssumx,nmax2,fnmax,"pair:dssumx");
864 memory->grow(dssumy,nmax2,fnmax,"pair:dssumy");
865 memory->grow(dssumz,nmax2,fnmax,"pair:dssumz");
866 }
867 }
868 for (j=0;j<f;j++) {
869 features[j]=0;
870 }
871 for (j=0;j<f*jnum;j++) {
872 dfeaturesx[j]=dfeaturesy[j]=dfeaturesz[j]=0;
873 }
874 //screening is calculated once for all atoms if any fingerprint uses it.
875 if (dospin) {
876 for (j=0;j<f*jnum;j++) {
877 sx[j]=sy[j]=sz[j]=0;
878 }
879 }
880 if (doscreen) screening(ii,0,jnum-1);
881 if (allscreen) screen_neighbor_list(&jnum);
882 //do fingerprints for atom type
883 len = fingerprintperelement[itype];
884 for (j=0;j<len;j++) {
885 if (fingerprints[itype][j]->spin==false && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,ii,nn,xn,yn,zn,tn,jnum-1,jl);
886 else if (fingerprints[itype][j]->spin==false && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl);
887 else if (fingerprints[itype][j]->spin==true && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,ii,nn,xn,yn,zn,tn,jnum-1,jl);
888 else if (fingerprints[itype][j]->spin==true && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl);
889 }
890 itype = nelements;
891 //do fingerprints for type "all"
892 len = fingerprintperelement[itype];
893 for (j=0;j<len;j++) {
894 if (fingerprints[itype][j]->spin==false && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,ii,nn,xn,yn,zn,tn,jnum-1,jl);
895 else if (fingerprints[itype][j]->spin==false && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl);
896 else if (fingerprints[itype][j]->spin==true && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,ii,nn,xn,yn,zn,tn,jnum-1,jl);
897 else if (fingerprints[itype][j]->spin==true && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl);
898 }
899 //run fingerprints through network
900 if (dospin) {
901 propagateforwardspin(&energy,force,fm,ii,jnum);
902 } else {
903 propagateforward(&energy,force,ii,jnum);
904 }
905 }
906 if (vflag_fdotr) virial_fdotr_compute();
907 delete[] sims;
908 }
909
cull_neighbor_list(int * jnum,int i,int sn)910 void PairRANN::cull_neighbor_list(int* jnum,int i,int sn) {
911 int *jlist,j,count,jj,*type,jtype;
912 double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
913 double **x = sims[sn].x;
914 xtmp = x[i][0];
915 ytmp = x[i][1];
916 ztmp = x[i][2];
917 type = sims[sn].type;
918 jlist = sims[sn].firstneigh[i];
919 count = 0;
920 for (jj=0;jj<jnum[0];jj++) {
921 j = jlist[jj];
922 j &= NEIGHMASK;
923 jtype = map[type[j]];
924 delx = xtmp - x[j][0];
925 dely = ytmp - x[j][1];
926 delz = ztmp - x[j][2];
927 rsq = delx*delx + dely*dely + delz*delz;
928 if (rsq>cutmax*cutmax) {
929 continue;
930 }
931 xn[count]=delx;
932 yn[count]=dely;
933 zn[count]=delz;
934 tn[count]=jtype;
935 jl[count]=j;
936 count++;
937 }
938 jnum[0]=count+1;
939 }
940
screen_neighbor_list(int * jnum)941 void PairRANN::screen_neighbor_list(int *jnum) {
942 int jj,kk,count,count1;
943 count = 0;
944 for (jj=0;jj<jnum[0]-1;jj++) {
945 if (Bij[jj]) {
946 count1 = 0;
947 if (jj!=count) {
948 xn[count]=xn[jj];
949 yn[count]=yn[jj];
950 zn[count]=zn[jj];
951 tn[count]=tn[jj];
952 jl[count]=jl[jj];
953 Sik[count]=Sik[jj];
954 dSikx[count]=dSikx[jj];
955 dSiky[count]=dSiky[jj];
956 dSikz[count]=dSikz[jj];
957 }
958 for (kk=0;kk<jnum[0]-1;kk++) {
959 if (Bij[kk]) {
960 dSijkxc[count][count1] = dSijkx[jj*(jnum[0]-1)+kk];
961 dSijkyc[count][count1] = dSijky[jj*(jnum[0]-1)+kk];
962 dSijkzc[count][count1] = dSijkz[jj*(jnum[0]-1)+kk];
963 count1++;
964 }
965 }
966 count++;
967 }
968 }
969 jnum[0]=count+1;
970 for (jj=0;jj<count;jj++) {
971 Bij[jj]=true;
972 for (kk=0;kk<count;kk++) {
973 dSijkx[jj*count+kk] = dSijkxc[jj][kk];
974 dSijky[jj*count+kk] = dSijkyc[jj][kk];
975 dSijkz[jj*count+kk] = dSijkzc[jj][kk];
976 }
977 }
978 }
979
980
screening(int ii,int sid,int jnum)981 void PairRANN::screening(int ii,int sid,int jnum)
982 {
983 //see Baskes, Materials Chemistry and Physics 50 (1997) 152-1.58
984 int i,jj,kk,itype,jtype,ktype;
985 double Sijk,Cijk,Cn,Cd,Dij,Dik,Djk,C,dfc,dC;
986 PairRANN::Simulation *sim = &sims[sid];
987 double delx,dely,delz,rij,delx2,dely2,delz2,rik,delx3,dely3,delz3,rjk;
988 i = sim->ilist[ii];
989 itype = map[sim->type[i]];
990 for (int jj=0;jj<jnum;jj++) {
991 Sik[jj]=1;
992 Bij[jj]=true;
993 dSikx[jj]=0;
994 dSiky[jj]=0;
995 dSikz[jj]=0;
996 for (kk=0;kk<jnum;kk++) {
997 dSijkx[jj*jnum+kk]=0;
998 dSijky[jj*jnum+kk]=0;
999 dSijkz[jj*jnum+kk]=0;
1000 }
1001 }
1002 for (kk=0;kk<jnum;kk++) {//outer sum over k in accordance with source, some others reorder to outer sum over jj
1003 if (Bij[kk]==false) {continue;}
1004 ktype = tn[kk];
1005 delx2 = xn[kk];
1006 dely2 = yn[kk];
1007 delz2 = zn[kk];
1008 rik = delx2*delx2+dely2*dely2+delz2*delz2;
1009 if (rik>cutmax*cutmax) {
1010 Bij[kk]= false;
1011 continue;
1012 }
1013 for (jj=0;jj<jnum;jj++) {
1014 if (jj==kk) {continue;}
1015 if (Bij[jj]==false) {continue;}
1016 jtype = tn[jj];
1017 delx = xn[jj];
1018 dely = yn[jj];
1019 delz = zn[jj];
1020 rij = delx*delx+dely*dely+delz*delz;
1021 if (rij>cutmax*cutmax) {
1022 Bij[jj] = false;
1023 continue;
1024 }
1025 delx3 = delx2-delx;
1026 dely3 = dely2-dely;
1027 delz3 = delz2-delz;
1028 rjk = delx3*delx3+dely3*dely3+delz3*delz3;
1029 if (rik+rjk<=rij) {continue;}//bond angle > 90 degrees
1030 if (rik+rij<=rjk) {continue;}//bond angle > 90 degrees
1031 double Cmax = screening_max[itype*nelements*nelements+jtype*nelements+ktype];
1032 double Cmin = screening_min[itype*nelements*nelements+jtype*nelements+ktype];
1033 double temp1 = rij-rik+rjk;
1034 Cn = temp1*temp1-4*rij*rjk;
1035 //a few expanded computations provided for clarity:
1036 //Cn = (rij-rik+rjk)*(rij-rik+rjk)-4*rij*rjk;
1037 //Cd = (rij-rjk)*(rij-rjk)-rik*rik;
1038 //Cijk = 1+2*(rik*rij+rik*rjk-rik*rik)/(rik*rik-(rij-rjk)*(rij-rjk));
1039 temp1 = rij-rjk;
1040 Cd = temp1*temp1-rik*rik;
1041 Cijk = Cn/Cd;
1042 C = (Cijk-Cmin)/(Cmax-Cmin);
1043 if (C>=1) {continue;}
1044 else if (C<=0) {
1045 Bij[kk]=false;
1046 break;
1047 }
1048 dC = Cmax-Cmin;
1049 dC *= dC;
1050 dC *= dC;
1051 temp1 = 1-C;
1052 temp1 *= temp1;
1053 temp1 *= temp1;
1054 Sijk = 1-temp1;
1055 Sijk *= Sijk;
1056 Dij = 4*rik*(Cn+4*rjk*(rij+rik-rjk))/Cd/Cd;
1057 Dik = -4*(rij*Cn+rjk*Cn+8*rij*rik*rjk)/Cd/Cd;
1058 Djk = 4*rik*(Cn+4*rij*(rik-rij+rjk))/Cd/Cd;
1059 temp1 = Cijk-Cmax;
1060 double temp2 = temp1*temp1;
1061 dfc = 8*temp1*temp2/(temp2*temp2-dC);
1062 Sik[kk] *= Sijk;
1063 dSijkx[kk*jnum+jj] = dfc*(delx*Dij-delx3*Djk);
1064 dSikx[kk] += dfc*(delx2*Dik+delx3*Djk);
1065 dSijky[kk*jnum+jj] = dfc*(dely*Dij-dely3*Djk);
1066 dSiky[kk] += dfc*(dely2*Dik+dely3*Djk);
1067 dSijkz[kk*jnum+jj] = dfc*(delz*Dij-delz3*Djk);
1068 dSikz[kk] += dfc*(delz2*Dik+delz3*Djk);
1069 }
1070 }
1071 }
1072
1073
1074 //Called by getproperties. Propagate features and dfeatures through network. Updates force and energy
propagateforward(double * energy,double ** force,int ii,int jnum)1075 void PairRANN::propagateforward(double *energy,double **force,int ii,int jnum) {
1076 int i,j,k,jj,j1,itype,i1;
1077 int *ilist;
1078 ilist = listfull->ilist;
1079 int *type = atom->type;
1080 i1=ilist[ii];
1081 itype = map[type[i1]];
1082 NNarchitecture net1 = net[itype];
1083 int L = net1.layers-1;
1084 //energy output with forces from analytical derivatives
1085 double dsum1;
1086 int f = net1.dimensions[0];
1087 for (i=0;i<net1.layers-1;i++) {
1088 for (j=0;j<net1.dimensions[i+1];j++) {
1089 //energy forward propagation
1090 sum[j]=0;
1091 for (k=0;k<net1.dimensions[i];k++) {
1092 if (i==0&&j==0) {
1093 layer[k]=features[k];
1094 }
1095 sum[j] += net1.Weights[i][j*net1.dimensions[i]+k]*layer[k];
1096 }
1097 sum[j] += net1.Biases[i][j];
1098 dsum1 = activation[itype][i]->dactivation_function(sum[j]);
1099 sum[j] = activation[itype][i]->activation_function(sum[j]);
1100 if (i==L-1) {
1101 energy[j] = sum[j];
1102 if (eflag_atom) eatom[i1]=sum[j];
1103 if (eflag_global) eng_vdwl +=sum[j];
1104 }
1105 //force propagation
1106 for (jj=0;jj<jnum;jj++) {
1107 dlayersumx[jj][j]=0;
1108 dlayersumy[jj][j]=0;
1109 dlayersumz[jj][j]=0;
1110 for (k=0;k<net1.dimensions[i];k++) {
1111 if (i==0&&j==0) {
1112 dlayerx[jj][k]=dfeaturesx[jj*f+k];
1113 dlayery[jj][k]=dfeaturesy[jj*f+k];
1114 dlayerz[jj][k]=dfeaturesz[jj*f+k];
1115 }
1116 double w1 = net1.Weights[i][j*net1.dimensions[i]+k];
1117 dlayersumx[jj][j] += w1*dlayerx[jj][k];
1118 dlayersumy[jj][j] += w1*dlayery[jj][k];
1119 dlayersumz[jj][j] += w1*dlayerz[jj][k];
1120 }
1121 dlayersumx[jj][j] = dsum1*dlayersumx[jj][j];
1122 dlayersumy[jj][j] = dsum1*dlayersumy[jj][j];
1123 dlayersumz[jj][j] = dsum1*dlayersumz[jj][j];
1124 if (i==L-1 && jj < (jnum-1)) {
1125 j1 = jl[jj];
1126 force[j1][0]+=dlayersumx[jj][j];
1127 force[j1][1]+=dlayersumy[jj][j];
1128 force[j1][2]+=dlayersumz[jj][j];
1129 }
1130 }
1131 if (i==L-1) {
1132 j1 = ilist[ii];
1133 jj = jnum-1;
1134 force[j1][0]+=dlayersumx[jj][j];
1135 force[j1][1]+=dlayersumy[jj][j];
1136 force[j1][2]+=dlayersumz[jj][j];
1137 }
1138 }
1139 //update values for next iteration
1140 for (j=0;j<net1.dimensions[i+1];j++) {
1141 layer[j]=sum[j];
1142 for (jj=0;jj<jnum;jj++) {
1143 dlayerx[jj][j] = dlayersumx[jj][j];
1144 dlayery[jj][j] = dlayersumy[jj][j];
1145 dlayerz[jj][j] = dlayersumz[jj][j];
1146 }
1147 }
1148 }
1149 }
1150
1151 //Called by getproperties. Propagate features and dfeatures through network. Updates force and energy
propagateforwardspin(double * energy,double ** force,double ** fm,int ii,int jnum)1152 void PairRANN::propagateforwardspin(double * energy,double **force,double **fm,int ii,int jnum) {
1153 int i,j,k,jj,j1,itype,i1;
1154 int *ilist;
1155 ilist = listfull->ilist;
1156 int *type = atom->type;
1157 i1=ilist[ii];
1158 itype = map[type[i1]];
1159 NNarchitecture net1 = net[itype];
1160 int L = net1.layers-1;
1161 //energy output with forces from analytical derivatives
1162 double dsum1;
1163 int f = net1.dimensions[0];
1164 for (i=0;i<net1.layers-1;i++) {
1165 for (j=0;j<net1.dimensions[i+1];j++) {
1166 //energy forward propagation
1167 sum[j]=0;
1168 for (k=0;k<net1.dimensions[i];k++) {
1169 if (i==0&&j==0) {
1170 layer[k]=features[k];
1171 }
1172 sum[j] += net1.Weights[i][j*net1.dimensions[i]+k]*layer[k];
1173 }
1174 sum[j] += net1.Biases[i][j];
1175 dsum1 = activation[itype][i]->dactivation_function(sum[j]);
1176 sum[j] = activation[itype][i]->activation_function(sum[j]);
1177 if (i==L-1) {
1178 energy[j] = sum[j];
1179 if (eflag_atom) eatom[i1]=sum[j];
1180 if (eflag_global) eng_vdwl +=sum[j];
1181 }
1182 //force propagation
1183 for (jj=0;jj<jnum;jj++) {
1184 dlayersumx[jj][j]=0;
1185 dlayersumy[jj][j]=0;
1186 dlayersumz[jj][j]=0;
1187 dssumx[jj][j]=0;
1188 dssumy[jj][j]=0;
1189 dssumz[jj][j]=0;
1190 for (k=0;k<net1.dimensions[i];k++) {
1191 if (i==0&&j==0) {
1192 dlayerx[jj][k]=dfeaturesx[jj*f+k];
1193 dlayery[jj][k]=dfeaturesy[jj*f+k];
1194 dlayerz[jj][k]=dfeaturesz[jj*f+k];
1195 dsx[jj][k]=-sx[jj*f+k];
1196 dsy[jj][k]=-sy[jj*f+k];
1197 dsz[jj][k]=-sz[jj*f+k];
1198 }
1199 double w1 = net1.Weights[i][j*net1.dimensions[i]+k];
1200 dlayersumx[jj][j] += w1*dlayerx[jj][k];
1201 dlayersumy[jj][j] += w1*dlayery[jj][k];
1202 dlayersumz[jj][j] += w1*dlayerz[jj][k];
1203 dssumx[jj][j] += w1*dsx[jj][k];
1204 dssumy[jj][j] += w1*dsy[jj][k];
1205 dssumz[jj][j] += w1*dsz[jj][k];
1206 }
1207 dlayersumx[jj][j] = dsum1*dlayersumx[jj][j];
1208 dlayersumy[jj][j] = dsum1*dlayersumy[jj][j];
1209 dlayersumz[jj][j] = dsum1*dlayersumz[jj][j];
1210 dssumx[jj][j] *= dsum1;
1211 dssumy[jj][j] *= dsum1;
1212 dssumz[jj][j] *= dsum1;
1213 if (i==L-1 && jj < (jnum-1)) {
1214 j1 = jl[jj];
1215 force[j1][0]+=dlayersumx[jj][j];
1216 force[j1][1]+=dlayersumy[jj][j];
1217 force[j1][2]+=dlayersumz[jj][j];
1218 fm[j1][0]+=dssumx[jj][j];
1219 fm[j1][1]+=dssumy[jj][j];
1220 fm[j1][2]+=dssumz[jj][j];
1221 }
1222 }
1223 if (i==L-1) {
1224 j1 = ilist[ii];
1225 jj = jnum-1;
1226 force[j1][0]+=dlayersumx[jj][j];
1227 force[j1][1]+=dlayersumy[jj][j];
1228 force[j1][2]+=dlayersumz[jj][j];
1229 fm[j1][0]+=dssumx[jj][j];
1230 fm[j1][1]+=dssumy[jj][j];
1231 fm[j1][2]+=dssumz[jj][j];
1232 }
1233 }
1234 //update values for next iteration
1235 for (j=0;j<net1.dimensions[i+1];j++) {
1236 layer[j]=sum[j];
1237 for (jj=0;jj<jnum;jj++) {
1238 dlayerx[jj][j] = dlayersumx[jj][j];
1239 dlayery[jj][j] = dlayersumy[jj][j];
1240 dlayerz[jj][j] = dlayersumz[jj][j];
1241 dsx[jj][j] = dssumx[jj][j];
1242 dsy[jj][j] = dssumy[jj][j];
1243 dsz[jj][j] = dssumz[jj][j];
1244 }
1245 }
1246 }
1247 }
1248
init_list(int,NeighList * ptr)1249 void PairRANN::init_list(int /*which*/, NeighList *ptr)
1250 {
1251 listfull = ptr;
1252 }
1253
init_style()1254 void PairRANN::init_style()
1255 {
1256 int irequest_full = neighbor->request(this,instance_me);
1257 neighbor->requests[irequest_full]->id = 1;
1258 neighbor->requests[irequest_full]->half = 0;
1259 neighbor->requests[irequest_full]->full = 1;
1260 }
1261
1262
1263 /* ----------------------------------------------------------------------
1264 init for one type pair i,j and corresponding j,i
1265 ------------------------------------------------------------------------- */
1266
init_one(int,int)1267 double PairRANN::init_one(int /*i*/, int /*j*/)
1268 {
1269 return cutmax;
1270 }
1271
errorf(const char * file,int line,const char * message)1272 void PairRANN::errorf(const char *file, int line, const char * message) {
1273 error->one(file,line,message);
1274 }
1275
factorial(int n)1276 int PairRANN::factorial(int n) {
1277 return round(MathSpecial::factorial(n));
1278 }
1279
create_fingerprint(const char * style)1280 RANN::Fingerprint *PairRANN::create_fingerprint(const char *style)
1281 {
1282 if (strcmp(style,"radial")==0) {
1283 return new RANN::Fingerprint_radial(this);
1284 }
1285 else if (strcmp(style,"radialscreened")==0) {
1286 return new RANN::Fingerprint_radialscreened(this);
1287 }
1288 else if (strcmp(style,"radialscreenedspin")==0) {
1289 return new RANN::Fingerprint_radialscreenedspin(this);
1290 }
1291 else if (strcmp(style,"radialspin")==0) {
1292 return new RANN::Fingerprint_radialspin(this);
1293 }
1294 else if (strcmp(style,"bond")==0) {
1295 return new RANN::Fingerprint_bond(this);
1296 }
1297 else if (strcmp(style,"bondscreened")==0) {
1298 return new RANN::Fingerprint_bondscreened(this);
1299 }
1300 else if (strcmp(style,"bondscreenedspin")==0) {
1301 return new RANN::Fingerprint_bondscreenedspin(this);
1302 }
1303 else if (strcmp(style,"bondspin")==0) {
1304 return new RANN::Fingerprint_bondspin(this);
1305 }
1306 error->one(FLERR,"Unknown fingerprint style {}",style);
1307 return nullptr;
1308 }
1309
1310
create_activation(const char * style)1311 RANN::Activation *PairRANN::create_activation(const char *style)
1312 {
1313 if (strcmp(style,"linear")==0) {
1314 return new RANN::Activation_linear(this);
1315 }
1316 else if (strcmp(style,"sigI")==0) {
1317 return new RANN::Activation_sigI(this);
1318 }
1319 error->one(FLERR,"Unknown activation style {}",style);
1320 return nullptr;
1321 }
1322
1323