diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst index 7573a2031a..e43e61cd00 100644 --- a/doc/src/compute_sna_atom.rst +++ b/doc/src/compute_sna_atom.rst @@ -33,7 +33,7 @@ Syntax * R_1, R_2,... = list of cutoff radii, one for each type (distance units) * w_1, w_2,... = list of neighbor weights, one for each type * zero or more keyword/value pairs may be appended -* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* +* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* .. parsed-literal:: @@ -59,6 +59,9 @@ Syntax *bikflag* value = *0* or *1* (only implemented for compute snap) *0* = per-atom bispectrum descriptors are summed over atoms *1* = per-atom bispectrum descriptors are not summed over atoms + *switchinnerflag* values = *rinnerlist* *drinnerlist* + *rinnerlist* = *ntypes* values of rinner (distance units) + *drinnerlist* = *ntypes* values of drinner (distance units) Examples """""""" @@ -70,6 +73,7 @@ Examples compute vb all sna/atom 1.4 0.95 6 2.0 1.0 compute snap all snap 1.4 0.95 6 2.0 1.0 compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 chem 2 0 1 + compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 switchinnerflag 1.1 1.3 0.5 0.6 Description """"""""""" @@ -308,6 +312,26 @@ the resulting bispectrum rows are :math:`B_{i,k}` instead of just :math:`B_k`. In this case, the entries in the final column for these rows are set to zero. +The keyword *switchinnerflag* activates an additional radial switching +function similar to :math:`f_c(r)` above, but acting to switch off +smoothly contributions from neighbor atoms at short separation distances. +This is useful when SNAP is used in combination with a simple +repulsive potential. The keyword is followed by the *ntypes* +values for :math:`r_{inner}` and the *ntypes* +values for :math:`\Delta r_{inner}`. For a neighbor atom at +distance :math:`r`, its contribution is scaled by a multiplicative +factor :math:`f_{inner}(r)` defined as follows: + +.. math:: + + = & 0, r \leq r_{inner} \\ + f_{inner}(r) = & \frac{1}{2}(1 - \cos(\pi \frac{r-r_{inner}}{\Delta r_{inner}})), r_{inner} < r \leq r_{inner} + \Delta r_{inner} \\ + = & 1, r > r_{inner} + \Delta r_{inner} + +The values of :math:`r_{inner}` and :math:`\Delta r_{inner}` are +the arithmetic means of the values for the central atom of type I +and the neighbor atom of type J. + .. note:: If you have a bonded system, then the settings of :doc:`special_bonds diff --git a/doc/src/pair_snap.rst b/doc/src/pair_snap.rst index 1bc17fa8c8..350561fe4d 100644 --- a/doc/src/pair_snap.rst +++ b/doc/src/pair_snap.rst @@ -135,7 +135,8 @@ with #) anywhere. Each non-blank non-comment line must contain one keyword/value pair. The required keywords are *rcutfac* and *twojmax*\ . Optional keywords are *rfac0*, *rmin0*, *switchflag*, *bzeroflag*, *quadraticflag*, *chemflag*, -*bnormflag*, *wselfallflag*, *chunksize*, and *parallelthresh*\ . +*bnormflag*, *wselfallflag*, *switchinnerflag*, +*rinner*, *drinner*, *chunksize*, and *parallelthresh*\ . The default values for these keywords are @@ -147,6 +148,7 @@ The default values for these keywords are * *chemflag* = 0 * *bnormflag* = 0 * *wselfallflag* = 0 +* *switchinnerflag* = 0 * *chunksize* = 32768 * *parallelthresh* = 8192 @@ -189,6 +191,16 @@ corresponding *K*-vector of linear coefficients for element which must equal the number of unique elements appearing in the LAMMPS pair_coeff command, to avoid ambiguity in the number of coefficients. +The keyword *switchinnerflag* activates an additional switching function +that smoothly turns off contributions to the SNAP potential from neighbor +atoms at short separations. If *switchinnerflag* is set to 1 then +the additional keywords *rinner* and *drinner* must also be provided. +Each of these is followed by *nelements* values, where *nelements* +is the number of unique elements appearing in appearing in the LAMMPS +pair_coeff command. The element order should correspond to the order +in which elements first appear in the pair_coeff command reading from +left to right. + The keywords *chunksize* and *parallelthresh* are only applicable when using the pair style *snap* with the KOKKOS package on GPUs and are ignored otherwise. The *chunksize* keyword controls the number of atoms diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index ac4dceed65..4e05cd5182 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1928,6 +1928,7 @@ mbt MBytes mc McLachlan +mcmoves md mdf mdi @@ -2228,6 +2229,7 @@ Nelement Nelements nelems nemd +netapp netcdf netstat Nettleton diff --git a/potentials/InP_JCPA2020.mliap b/potentials/InP_JCPA2020.mliap index 19d044df97..f5d86590b6 100644 --- a/potentials/InP_JCPA2020.mliap +++ b/potentials/InP_JCPA2020.mliap @@ -1,4 +1,4 @@ -# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020) +# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, 124 5456 (2020) # Definition of SNAP+ZBL potential. diff --git a/potentials/InP_JCPA2020.mliap.descriptor b/potentials/InP_JCPA2020.mliap.descriptor index 287eac59c4..ed70da5403 100644 --- a/potentials/InP_JCPA2020.mliap.descriptor +++ b/potentials/InP_JCPA2020.mliap.descriptor @@ -1,4 +1,4 @@ -# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020) +# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, 124 5456 (2020) # required rcutfac 1.0 diff --git a/potentials/InP_JCPA2020.mliap.model b/potentials/InP_JCPA2020.mliap.model index 7ef9039ba2..9217ed4df7 100644 --- a/potentials/InP_JCPA2020.mliap.model +++ b/potentials/InP_JCPA2020.mliap.model @@ -1,4 +1,4 @@ -# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020) +# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, 124 5456 (2020) 2 241 0.000000000000 # B[0] Block = 1 Type = In diff --git a/potentials/InP_JCPA2020.snap b/potentials/InP_JCPA2020.snap index 1af0008b6f..90ff5d37ff 100644 --- a/potentials/InP_JCPA2020.snap +++ b/potentials/InP_JCPA2020.snap @@ -1,4 +1,4 @@ -# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020) +# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, 124 5456 (2020) # Definition of SNAP+ZBL potential. diff --git a/potentials/InP_JCPA2020.snapcoeff b/potentials/InP_JCPA2020.snapcoeff index d4be08e4e2..8a07006480 100644 --- a/potentials/InP_JCPA2020.snapcoeff +++ b/potentials/InP_JCPA2020.snapcoeff @@ -1,4 +1,4 @@ -# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020) +# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, 124 5456 (2020) 2 241 In 3.81205 1 diff --git a/potentials/InP_JCPA2020.snapparam b/potentials/InP_JCPA2020.snapparam index 0e764ac7ca..b699748032 100644 --- a/potentials/InP_JCPA2020.snapparam +++ b/potentials/InP_JCPA2020.snapparam @@ -1,4 +1,4 @@ -# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020) +# DATE: 2020-06-01 UNITS: metal CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, 124 5456 (2020) # required rcutfac 1.0 diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index 5475117a99..dd9bf5d812 100644 --- a/src/ML-IAP/mliap_descriptor_snap.cpp +++ b/src/ML-IAP/mliap_descriptor_snap.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -25,6 +24,7 @@ #include "mliap_data.h" #include "pair_mliap.h" #include "sna.h" +#include "tokenizer.h" #include #include @@ -36,17 +36,18 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename): - MLIAPDescriptor(lmp) +MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *_lmp, char *paramfilename) : MLIAPDescriptor(_lmp) { radelem = nullptr; wjelem = nullptr; snaptr = nullptr; + rinnerelem = nullptr; + drinnerelem = nullptr; + read_paramfile(paramfilename); - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); ndescriptors = snaptr->ncoeff; } @@ -58,13 +59,18 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP() memory->destroy(radelem); memory->destroy(wjelem); delete snaptr; + + if (switchinnerflag) { + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); + } } /* ---------------------------------------------------------------------- compute descriptors for each atom ---------------------------------------------------------------------- */ -void MLIAPDescriptorSNAP::compute_descriptors(class MLIAPData* data) +void MLIAPDescriptorSNAP::compute_descriptors(class MLIAPData *data) { int ij = 0; for (int ii = 0; ii < data->nlistatoms; ii++) { @@ -87,7 +93,11 @@ void MLIAPDescriptorSNAP::compute_descriptors(class MLIAPData* data) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); - snaptr->element[ninside] = jelem; // element index for chem snap + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5 * (rinnerelem[ielem] + rinnerelem[jelem]); + snaptr->drinnerij[ninside] = 0.5 * (drinnerelem[ielem] + drinnerelem[jelem]); + } + if (chemflag) snaptr->element[ninside] = jelem; ninside++; ij++; } @@ -106,14 +116,13 @@ void MLIAPDescriptorSNAP::compute_descriptors(class MLIAPData* data) for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) data->descriptors[ii][icoeff] = snaptr->blist[icoeff]; } - } /* ---------------------------------------------------------------------- compute forces for each atom ---------------------------------------------------------------------- */ -void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data) +void MLIAPDescriptorSNAP::compute_forces(class MLIAPData *data) { double fij[3]; double **f = atom->f; @@ -140,7 +149,11 @@ void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); - snaptr->element[ninside] = jelem; // element index for chem snap + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5 * (rinnerelem[ielem] + rinnerelem[jelem]); + snaptr->drinnerij[ninside] = 0.5 * (drinnerelem[ielem] + drinnerelem[jelem]); + } + if (chemflag) snaptr->element[ninside] = jelem; ninside++; ij++; } @@ -160,12 +173,7 @@ void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data) for (int jj = 0; jj < ninside; jj++) { int j = snaptr->inside[jj]; - if (chemflag) - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, snaptr->element[jj]); - else - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, 0); + snaptr->compute_duidrj(jj); snaptr->compute_deidrj(fij); @@ -179,18 +187,16 @@ void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data) // add in global and per-atom virial contributions // this is optional and has no effect on force calculation - if (data->vflag) - data->pairmliap->v_tally(i,j,fij,snaptr->rij[jj]); + if (data->vflag) data->pairmliap->v_tally(i, j, fij, snaptr->rij[jj]); } } - } /* ---------------------------------------------------------------------- calculate gradients of forces w.r.t. parameters ---------------------------------------------------------------------- */ -void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data) +void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData *data) { int ij = 0; for (int ii = 0; ii < data->nlistatoms; ii++) { @@ -215,7 +221,11 @@ void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); - snaptr->element[ninside] = jelem; // element index for chem snap + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5 * (rinnerelem[ielem] + rinnerelem[jelem]); + snaptr->drinnerij[ninside] = 0.5 * (drinnerelem[ielem] + drinnerelem[jelem]); + } + if (chemflag) snaptr->element[ninside] = jelem; ninside++; ij++; } @@ -234,13 +244,7 @@ void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data) for (int jj = 0; jj < ninside; jj++) { const int j = snaptr->inside[jj]; - if (chemflag) - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, snaptr->element[jj]); - else - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, 0); - + snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate gamma_lk*dB_k/dRi, -gamma_lk**dB_k/dRj @@ -248,24 +252,22 @@ void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data) for (int inz = 0; inz < data->gamma_nnz; inz++) { const int l = data->gamma_row_index[ii][inz]; const int k = data->gamma_col_index[ii][inz]; - data->gradforce[i][l] += data->gamma[ii][inz]*snaptr->dblist[k][0]; - data->gradforce[i][l+data->yoffset] += data->gamma[ii][inz]*snaptr->dblist[k][1]; - data->gradforce[i][l+data->zoffset] += data->gamma[ii][inz]*snaptr->dblist[k][2]; - data->gradforce[j][l] -= data->gamma[ii][inz]*snaptr->dblist[k][0]; - data->gradforce[j][l+data->yoffset] -= data->gamma[ii][inz]*snaptr->dblist[k][1]; - data->gradforce[j][l+data->zoffset] -= data->gamma[ii][inz]*snaptr->dblist[k][2]; + data->gradforce[i][l] += data->gamma[ii][inz] * snaptr->dblist[k][0]; + data->gradforce[i][l + data->yoffset] += data->gamma[ii][inz] * snaptr->dblist[k][1]; + data->gradforce[i][l + data->zoffset] += data->gamma[ii][inz] * snaptr->dblist[k][2]; + data->gradforce[j][l] -= data->gamma[ii][inz] * snaptr->dblist[k][0]; + data->gradforce[j][l + data->yoffset] -= data->gamma[ii][inz] * snaptr->dblist[k][1]; + data->gradforce[j][l + data->zoffset] -= data->gamma[ii][inz] * snaptr->dblist[k][2]; } - } } - } /* ---------------------------------------------------------------------- compute descriptor gradients for each neighbor atom ---------------------------------------------------------------------- */ -void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data) +void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData *data) { int ij = 0; for (int ii = 0; ii < data->nlistatoms; ii++) { @@ -289,7 +291,11 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); - snaptr->element[ninside] = jelem; // element index for chem snap + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5 * (rinnerelem[ielem] + rinnerelem[jelem]); + snaptr->drinnerij[ninside] = 0.5 * (drinnerelem[ielem] + drinnerelem[jelem]); + } + if (chemflag) snaptr->element[ninside] = jelem; ninside++; ij++; } @@ -307,13 +313,7 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data) ij = ij0; for (int jj = 0; jj < ninside; jj++) { - if (chemflag) - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, snaptr->element[jj]); - else - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, 0); - + snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate dB_k^i/dRi, dB_k^i/dRj @@ -326,7 +326,6 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data) ij++; } } - } /* ---------------------------------------------------------------------- @@ -361,6 +360,12 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 0; + + // set local input checks + + int rinnerflag = 0; + int drinnerflag = 0; for (int i = 0; i < nelements; i++) delete[] elements[i]; delete[] elements; @@ -372,128 +377,138 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) FILE *fpparam; if (comm->me == 0) { - fpparam = utils::open_potential(paramfilename,lmp,nullptr); + fpparam = utils::open_potential(paramfilename, lmp, nullptr); if (fpparam == nullptr) - error->one(FLERR,"Cannot open SNAP parameter file {}: {}", - paramfilename, utils::getsyserror()); + error->one(FLERR, "Cannot open SNAP parameter file {}: {}", paramfilename, + utils::getsyserror()); } - char line[MAXLINE],*ptr; + char line[MAXLINE], *ptr; int eof = 0; - int n,nwords; + int n; while (true) { if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fpparam); + ptr = fgets(line, MAXLINE, fpparam); if (ptr == nullptr) { eof = 1; fclose(fpparam); - } else n = strlen(line) + 1; + } else + n = strlen(line) + 1; } - MPI_Bcast(&eof,1,MPI_INT,0,world); + MPI_Bcast(&eof, 1, MPI_INT, 0, world); if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); + MPI_Bcast(&n, 1, MPI_INT, 0, world); + MPI_Bcast(line, n, MPI_CHAR, 0, world); // strip comment, skip line if blank - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - nwords = utils::count_words(line); - if (nwords == 0) continue; + if ((ptr = strchr(line, '#'))) *ptr = '\0'; - // words = ptrs to all words in line // strip single and double quotes from words + Tokenizer words(line, "\"' \t\n\t\f"); + if (words.count() == 0) continue; - char* keywd = strtok(line,"' \t\n\r\f"); - char* keyval = strtok(nullptr,"' \t\n\r\f"); - - if (comm->me == 0) - utils::logmesg(lmp,"SNAP keyword {} {} \n", keywd, keyval); + auto keywd = words.next(); // check for keywords with one value per element - if (strcmp(keywd,"elems") == 0 || - strcmp(keywd,"radelems") == 0 || - strcmp(keywd,"welems") == 0) { + if ((keywd == "elems") || (keywd == "radelems") || (keywd == "welems") || + (keywd == "rinnerelems") || (keywd == "drinnerelems")) { - if (nelementsflag == 0 || nwords != nelements+1) - error->all(FLERR,"Incorrect SNAP parameter file"); + if ((nelementsflag == 0) || (words.count() != nelements + 1)) + error->all(FLERR, "Incorrect SNAP parameter file"); - if (strcmp(keywd,"elems") == 0) { - for (int ielem = 0; ielem < nelements; ielem++) { - elements[ielem] = utils::strdup(keyval); - keyval = strtok(nullptr,"' \t\n\r\f"); - } + if (comm->me == 0) utils::logmesg(lmp, "SNAP keyword {} \n", utils::trim(line)); + + if (keywd == "elems") { + for (int ielem = 0; ielem < nelements; ielem++) + elements[ielem] = utils::strdup(words.next()); elementsflag = 1; - } else if (strcmp(keywd,"radelems") == 0) { - for (int ielem = 0; ielem < nelements; ielem++) { - radelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); - keyval = strtok(nullptr,"' \t\n\r\f"); - } + } else if (keywd == "radelems") { + for (int ielem = 0; ielem < nelements; ielem++) + radelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); radelemflag = 1; - } else if (strcmp(keywd,"welems") == 0) { - for (int ielem = 0; ielem < nelements; ielem++) { - wjelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); - keyval = strtok(nullptr,"' \t\n\r\f"); - } + } else if (keywd == "welems") { + for (int ielem = 0; ielem < nelements; ielem++) + wjelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); wjelemflag = 1; + } else if (keywd == "rinnerelems") { + for (int ielem = 0; ielem < nelements; ielem++) + rinnerelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); + rinnerflag = 1; + } else if (keywd == "drinnerelems") { + for (int ielem = 0; ielem < nelements; ielem++) + drinnerelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); + rinnerflag = 1; } } else { - // all other keywords take one value + // all other keywords take one value - if (nwords != 2) - error->all(FLERR,"Incorrect SNAP parameter file"); + if (words.count() != 2) error->all(FLERR, "Incorrect SNAP parameter file"); - if (strcmp(keywd,"nelems") == 0) { - nelements = atoi(keyval); - elements = new char*[nelements]; - memory->create(radelem,nelements,"mliap_snap_descriptor:radelem"); - memory->create(wjelem,nelements,"mliap_snap_descriptor:wjelem"); + auto keyval = words.next(); + if (comm->me == 0) utils::logmesg(lmp, "SNAP keyword {} {} \n", keywd, keyval); + + if (keywd == "nelems") { + nelements = utils::inumeric(FLERR, keyval, false, lmp); + elements = new char *[nelements]; + memory->create(radelem, nelements, "mliap_snap_descriptor:radelem"); + memory->create(wjelem, nelements, "mliap_snap_descriptor:wjelem"); + memory->create(rinnerelem, nelements, "mliap_snap_descriptor:rinner"); + memory->create(drinnerelem, nelements, "mliap_snap_descriptor:drinner"); nelementsflag = 1; - } else if (strcmp(keywd,"rcutfac") == 0) { - rcutfac = atof(keyval); + } else if (keywd == "rcutfac") { + rcutfac = utils::numeric(FLERR, keyval, false, lmp); rcutfacflag = 1; - } else if (strcmp(keywd,"twojmax") == 0) { - twojmax = atoi(keyval); + } else if (keywd == "twojmax") { + twojmax = utils::inumeric(FLERR, keyval, false, lmp); twojmaxflag = 1; - } else if (strcmp(keywd,"rfac0") == 0) - rfac0 = atof(keyval); - else if (strcmp(keywd,"rmin0") == 0) - rmin0 = atof(keyval); - else if (strcmp(keywd,"switchflag") == 0) - switchflag = atoi(keyval); - else if (strcmp(keywd,"bzeroflag") == 0) - bzeroflag = atoi(keyval); - else if (strcmp(keywd,"chemflag") == 0) - chemflag = atoi(keyval); - else if (strcmp(keywd,"bnormflag") == 0) - bnormflag = atoi(keyval); - else if (strcmp(keywd,"wselfallflag") == 0) - wselfallflag = atoi(keyval); + } else if (keywd == "rfac0") + rfac0 = utils::numeric(FLERR, keyval, false, lmp); + else if (keywd == "rmin0") + rmin0 = utils::numeric(FLERR, keyval, false, lmp); + else if (keywd == "switchflag") + switchflag = utils::inumeric(FLERR, keyval, false, lmp); + else if (keywd == "bzeroflag") + bzeroflag = utils::inumeric(FLERR, keyval, false, lmp); + else if (keywd == "chemflag") + chemflag = utils::inumeric(FLERR, keyval, false, lmp); + else if (keywd == "bnormflag") + bnormflag = utils::inumeric(FLERR, keyval, false, lmp); + else if (keywd == "wselfallflag") + wselfallflag = utils::inumeric(FLERR, keyval, false, lmp); + else if (keywd == "switchinnerflag") + switchinnerflag = utils::inumeric(FLERR, keyval, false, lmp); else - error->all(FLERR,"Incorrect SNAP parameter file"); - + error->all(FLERR, "Incorrect SNAP parameter file"); } } - if (!rcutfacflag || !twojmaxflag || !nelementsflag || - !elementsflag || !radelemflag || !wjelemflag) - error->all(FLERR,"Incorrect SNAP parameter file"); + if (!rcutfacflag || !twojmaxflag || !nelementsflag || !elementsflag || !radelemflag || + !wjelemflag) + error->all(FLERR, "Incorrect SNAP parameter file"); + + if (switchinnerflag && !(rinnerflag && drinnerflag)) + error->all(FLERR, "Incorrect SNAP parameter file"); + + if (!switchinnerflag && (rinnerflag || drinnerflag)) + error->all(FLERR, "Incorrect SNAP parameter file"); // construct cutsq double cut; cutmax = 0.0; - memory->create(cutsq,nelements,nelements,"mliap/descriptor/snap:cutsq"); + memory->create(cutsq, nelements, nelements, "mliap/descriptor/snap:cutsq"); for (int ielem = 0; ielem < nelements; ielem++) { - cut = 2.0*radelem[ielem]*rcutfac; + cut = 2.0 * radelem[ielem] * rcutfac; if (cut > cutmax) cutmax = cut; - cutsq[ielem][ielem] = cut*cut; - for (int jelem = ielem+1; jelem < nelements; jelem++) { - cut = (radelem[ielem]+radelem[jelem])*rcutfac; - cutsq[ielem][jelem] = cutsq[jelem][ielem] = cut*cut; + cutsq[ielem][ielem] = cut * cut; + for (int jelem = ielem + 1; jelem < nelements; jelem++) { + cut = (radelem[ielem] + radelem[jelem]) * rcutfac; + cutsq[ielem][jelem] = cutsq[jelem][ielem] = cut * cut; } } } @@ -505,7 +520,7 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) double MLIAPDescriptorSNAP::memory_usage() { double bytes = MLIAPDescriptor::memory_usage(); - bytes += snaptr->memory_usage(); // SNA object + bytes += snaptr->memory_usage(); // SNA object return bytes; } diff --git a/src/ML-IAP/mliap_descriptor_snap.h b/src/ML-IAP/mliap_descriptor_snap.h index 9098015a4f..2136eefc5e 100644 --- a/src/ML-IAP/mliap_descriptor_snap.h +++ b/src/ML-IAP/mliap_descriptor_snap.h @@ -39,7 +39,11 @@ class MLIAPDescriptorSNAP : public MLIAPDescriptor { int twojmax, switchflag, bzeroflag; int chemflag, bnormflag, wselfallflag; + int switchinnerflag; double rfac0, rmin0; + + double* rinnerelem; + double* drinnerelem; }; } // namespace LAMMPS_NS diff --git a/src/ML-IAP/mliap_model_nn.cpp b/src/ML-IAP/mliap_model_nn.cpp index 24fddda766..a68b29eefa 100644 --- a/src/ML-IAP/mliap_model_nn.cpp +++ b/src/ML-IAP/mliap_model_nn.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -33,8 +32,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -MLIAPModelNN::MLIAPModelNN(LAMMPS* lmp, char* coefffilename) : - MLIAPModel(lmp, coefffilename) +MLIAPModelNN::MLIAPModelNN(LAMMPS *_lmp, char *coefffilename) : MLIAPModel(_lmp, coefffilename) { nnodes = nullptr; activation = nullptr; @@ -47,9 +45,9 @@ MLIAPModelNN::MLIAPModelNN(LAMMPS* lmp, char* coefffilename) : MLIAPModelNN::~MLIAPModelNN() { - memory->destroy(nnodes); - memory->destroy(activation); - memory->destroy(scale); + memory->destroy(nnodes); + memory->destroy(activation); + memory->destroy(scale); } /* ---------------------------------------------------------------------- @@ -59,7 +57,7 @@ MLIAPModelNN::~MLIAPModelNN() int MLIAPModelNN::get_nparams() { if (nparams == 0) - if (ndescriptors == 0) error->all(FLERR,"ndescriptors not defined"); + if (ndescriptors == 0) error->all(FLERR, "ndescriptors not defined"); return nparams; } @@ -71,10 +69,10 @@ void MLIAPModelNN::read_coeffs(char *coefffilename) FILE *fpcoeff; if (comm->me == 0) { - fpcoeff = utils::open_potential(coefffilename,lmp,nullptr); + fpcoeff = utils::open_potential(coefffilename, lmp, nullptr); if (fpcoeff == nullptr) - error->one(FLERR,"Cannot open MLIAPModel coeff file {}: {}", - coefffilename,utils::getsyserror()); + error->one(FLERR, "Cannot open MLIAPModel coeff file {}: {}", coefffilename, + utils::getsyserror()); } char line[MAXLINE], *ptr, *tstr; @@ -84,24 +82,24 @@ void MLIAPModelNN::read_coeffs(char *coefffilename) int nwords = 0; while (nwords == 0) { if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fpcoeff); + ptr = fgets(line, MAXLINE, fpcoeff); if (ptr == nullptr) { eof = 1; fclose(fpcoeff); - } else n = strlen(line) + 1; + } else + n = strlen(line) + 1; } - MPI_Bcast(&eof,1,MPI_INT,0,world); + MPI_Bcast(&eof, 1, MPI_INT, 0, world); if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); + MPI_Bcast(&n, 1, MPI_INT, 0, world); + MPI_Bcast(line, n, MPI_CHAR, 0, world); // strip comment, skip line if blank - if ((ptr = strchr(line,'#'))) *ptr = '\0'; + if ((ptr = strchr(line, '#'))) *ptr = '\0'; nwords = utils::count_words(line); } - if (nwords != 2) - error->all(FLERR,"Incorrect format in MLIAPModel coefficient file"); + if (nwords != 2) error->all(FLERR, "Incorrect format in MLIAPModel coefficient file"); // words = ptrs to all words in line // strip single and double quotes from words @@ -111,14 +109,13 @@ void MLIAPModelNN::read_coeffs(char *coefffilename) nelements = coeffs.next_int(); nparams = coeffs.next_int(); } catch (TokenizerException &e) { - error->all(FLERR,"Incorrect format in MLIAPModel coefficient " - "file: {}",e.what()); + error->all(FLERR, "Incorrect format in MLIAPModel coefficient file: {}", e.what()); } // set up coeff lists memory->destroy(coeffelem); - memory->create(coeffelem,nelements,nparams,"mliap_snap_model:coeffelem"); + memory->create(coeffelem, nelements, nparams, "mliap_snap_model:coeffelem"); int stats = 0; int ielem = 0; @@ -126,85 +123,94 @@ void MLIAPModelNN::read_coeffs(char *coefffilename) while (true) { if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fpcoeff); + ptr = fgets(line, MAXLINE, fpcoeff); if (ptr == nullptr) { eof = 1; fclose(fpcoeff); - } else n = strlen(line) + 1; + } else + n = strlen(line) + 1; } - MPI_Bcast(&eof,1,MPI_INT,0,world); + MPI_Bcast(&eof, 1, MPI_INT, 0, world); if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); + MPI_Bcast(&n, 1, MPI_INT, 0, world); + MPI_Bcast(line, n, MPI_CHAR, 0, world); // strip comment, skip line if blank - if ((ptr = strchr(line,'#'))) *ptr = '\0'; + if ((ptr = strchr(line, '#'))) *ptr = '\0'; + nwords = utils::trim_and_count_words(line); if (nwords == 0) continue; - if (stats == 0) { // Header NET - tstr = strtok(line,"' \t\n\r\f"); - if (strncmp(tstr, "NET", 3) != 0) error->all(FLERR,"Incorrect format in MLIAPModel coefficient file"); + ValueTokenizer values(line, "\"' \t\n\t\f"); + if (stats == 0) { // Header NET + auto tstr = values.next_string(); + if (tstr.substr(0, 3) != "NET") + error->all(FLERR, "Incorrect format in MLIAPModel coefficient file"); - ndescriptors = atoi(strtok(nullptr,"' \t\n\r\f")); - nlayers = atoi(strtok(nullptr,"' \t\n\r\f")); + ndescriptors = values.next_int(); + nlayers = values.next_int(); - memory->create(activation,nlayers,"mliap_model:activation"); - memory->create(nnodes,nlayers,"mliap_model:nnodes"); - memory->create(scale,nelements,2,ndescriptors,"mliap_model:scale"); + memory->create(activation, nlayers, "mliap_model:activation"); + memory->create(nnodes, nlayers, "mliap_model:nnodes"); + memory->create(scale, nelements, 2, ndescriptors, "mliap_model:scale"); for (int ilayer = 0; ilayer < nlayers; ilayer++) { - tstr = strtok(nullptr,"' \t\n\r\f"); - nnodes[ilayer] = atoi(strtok(nullptr,"' \t\n\r\f")); + tstr = values.next_string(); + nnodes[ilayer] = values.next_int(); - if (strncmp(tstr, "linear", 6) == 0) activation[ilayer] = 0; - else if (strncmp(tstr, "sigmoid", 7) == 0) activation[ilayer] = 1; - else if (strncmp(tstr, "tanh", 4) == 0) activation[ilayer] = 2; - else if (strncmp(tstr, "relu", 4) == 0) activation[ilayer] = 3; - else activation[ilayer] = 4; + if (tstr == "linear") + activation[ilayer] = 0; + else if (tstr == "sigmoid") + activation[ilayer] = 1; + else if (tstr == "tanh") + activation[ilayer] = 2; + else if (tstr == "relu") + activation[ilayer] = 3; + else + activation[ilayer] = 4; } stats = 1; } else if (stats == 1) { - scale[ielem][0][l] = atof(strtok(line,"' \t\n\r\f")); - for (int icoeff = 1; icoeff < nwords; icoeff++) { - scale[ielem][0][l+icoeff] = atof(strtok(nullptr,"' \t\n\r\f")); - } - l += nwords; - if (l == ndescriptors) { - stats = 2; - l = 0; - } + scale[ielem][0][l] = values.next_double(); + for (int icoeff = 1; icoeff < nwords; icoeff++) { + scale[ielem][0][l + icoeff] = values.next_double(); + } + l += nwords; + if (l == ndescriptors) { + stats = 2; + l = 0; + } } else if (stats == 2) { - scale[ielem][1][l] = atof(strtok(line,"' \t\n\r\f")); - for (int icoeff = 1; icoeff < nwords; icoeff++) { - scale[ielem][1][l+icoeff] = atof(strtok(nullptr,"' \t\n\r\f")); - } - l += nwords; - if (l == ndescriptors) { - stats = 3; - l = 0; - } + scale[ielem][1][l] = values.next_double(); + for (int icoeff = 1; icoeff < nwords; icoeff++) { + scale[ielem][1][l + icoeff] = values.next_double(); + } + l += nwords; + if (l == ndescriptors) { + stats = 3; + l = 0; + } - // set up coeff lists + // set up coeff lists } else if (stats == 3) { - if (nwords > 30) error->all(FLERR,"Incorrect format in MLIAPModel coefficient file"); + if (nwords > 30) error->all(FLERR, "Incorrect format in MLIAPModel coefficient file"); - coeffelem[ielem][l] = atof(strtok(line,"' \t\n\r\f")); - for (int icoeff = 1; icoeff < nwords; icoeff++) { - coeffelem[ielem][l+icoeff] = atof(strtok(nullptr,"' \t\n\r\f")); - } - l += nwords; - if (l == nparams) { - stats = 1; - l = 0; - ielem++; - } + coeffelem[ielem][l] = values.next_double(); + for (int icoeff = 1; icoeff < nwords; icoeff++) { + coeffelem[ielem][l + icoeff] = values.next_double(); + } + l += nwords; + if (l == nparams) { + stats = 1; + l = 0; + ielem++; + } } } } @@ -214,153 +220,127 @@ void MLIAPModelNN::read_coeffs(char *coefffilename) for each atom beta_i = dE(B_i)/dB_i ---------------------------------------------------------------------- */ -void MLIAPModelNN::compute_gradients(MLIAPData* data) +void MLIAPModelNN::compute_gradients(MLIAPData *data) { data->energy = 0.0; for (int ii = 0; ii < data->nlistatoms; ii++) { - const int ielem = data->ielems[ii]; - const int nl = nlayers; + const int ielem = data->ielems[ii]; + const int nl = nlayers; - double* coeffi = coeffelem[ielem]; - double** scalei = scale[ielem]; - double **nodes, **dnodes, **bnodes; + double *coeffi = coeffelem[ielem]; + double **scalei = scale[ielem]; + double **nodes, **dnodes, **bnodes; - nodes = new double*[nl]; - dnodes = new double*[nl]; - bnodes = new double*[nl]; - for (int l=0; lndescriptors; icoeff++) { - nodes[0][n] += coeffi[n*((data->ndescriptors)+1)+icoeff+1] * - (data->descriptors[ii][icoeff] - scalei[0][icoeff]) / - scalei[1][icoeff]; - } - if (activation[0] == 1) { - nodes[0][n] = sigm(nodes[0][n] + - coeffi[n*((data->ndescriptors)+1)], - dnodes[0][n]); - } - else if (activation[0] == 2) { - nodes[0][n] = tanh(nodes[0][n] + - coeffi[n*((data->ndescriptors)+1)], - dnodes[0][n]); - } - else if (activation[0] == 3) { - nodes[0][n] = relu(nodes[0][n] + - coeffi[n*((data->ndescriptors)+1)], - dnodes[0][n]); - } - else { - nodes[0][n] += coeffi[n*((data->ndescriptors)+1)]; - dnodes[0][n] = 1; - } - } - - // hidden~output - - int k = 0; - if (nl > 1) { - k += ((data->ndescriptors)+1)*nnodes[0]; - for (int l=1; l < nl; l++) { - for (int n=0; n < nnodes[l]; n++) { - nodes[l][n] = 0; - for (int j=0; j < nnodes[l-1]; j++) { - nodes[l][n] += coeffi[k+n*(nnodes[l-1]+1)+j+1] * - nodes[l-1][j]; - } - if (activation[l] == 1) { - nodes[l][n] = sigm(nodes[l][n] + - coeffi[k+n*(nnodes[l-1]+1)], - dnodes[l][n]); - } - else if (activation[l] == 2) { - nodes[l][n] = tanh(nodes[l][n] + - coeffi[k+n*(nnodes[l-1]+1)], - dnodes[l][n]); - } - else if (activation[l] == 3) { - nodes[l][n] = relu(nodes[l][n] + - coeffi[k+n*(nnodes[l-1]+1)], - dnodes[l][n]); - } - else { - nodes[l][n] += coeffi[k+n*(nnodes[l-1]+1)]; - dnodes[l][n] = 1; - } - } - k += (nnodes[l-1]+1)*nnodes[l]; - } - } - - // backwardprop - // output layer dnode initialized to 1. - - for (int n=0; n 1) { - for (int l=nl-1; l>0; l--) { - k -= (nnodes[l-1]+1)*nnodes[l]; - for (int n=0; n= 1) { - bnodes[l-1][n] *= dnodes[l-1][n]; - } - } - } - } + // forwardprop + // input - hidden1 + for (int n = 0; n < nnodes[0]; n++) { + nodes[0][n] = 0; for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { - data->betas[ii][icoeff] = 0; - for (int j=0; jbetas[ii][icoeff] += - coeffi[j*((data->ndescriptors)+1)+icoeff+1] * - bnodes[0][j]; - } - data->betas[ii][icoeff] = data->betas[ii][icoeff]/scalei[1][icoeff]; + nodes[0][n] += coeffi[n * ((data->ndescriptors) + 1) + icoeff + 1] * + (data->descriptors[ii][icoeff] - scalei[0][icoeff]) / scalei[1][icoeff]; } + if (activation[0] == 1) { + nodes[0][n] = sigm(nodes[0][n] + coeffi[n * ((data->ndescriptors) + 1)], dnodes[0][n]); + } else if (activation[0] == 2) { + nodes[0][n] = tanh(nodes[0][n] + coeffi[n * ((data->ndescriptors) + 1)], dnodes[0][n]); + } else if (activation[0] == 3) { + nodes[0][n] = relu(nodes[0][n] + coeffi[n * ((data->ndescriptors) + 1)], dnodes[0][n]); + } else { + nodes[0][n] += coeffi[n * ((data->ndescriptors) + 1)]; + dnodes[0][n] = 1; + } + } - if (data->eflag) { + // hidden~output + + int k = 0; + if (nl > 1) { + k += ((data->ndescriptors) + 1) * nnodes[0]; + for (int l = 1; l < nl; l++) { + for (int n = 0; n < nnodes[l]; n++) { + nodes[l][n] = 0; + for (int j = 0; j < nnodes[l - 1]; j++) { + nodes[l][n] += coeffi[k + n * (nnodes[l - 1] + 1) + j + 1] * nodes[l - 1][j]; + } + if (activation[l] == 1) { + nodes[l][n] = sigm(nodes[l][n] + coeffi[k + n * (nnodes[l - 1] + 1)], dnodes[l][n]); + } else if (activation[l] == 2) { + nodes[l][n] = tanh(nodes[l][n] + coeffi[k + n * (nnodes[l - 1] + 1)], dnodes[l][n]); + } else if (activation[l] == 3) { + nodes[l][n] = relu(nodes[l][n] + coeffi[k + n * (nnodes[l - 1] + 1)], dnodes[l][n]); + } else { + nodes[l][n] += coeffi[k + n * (nnodes[l - 1] + 1)]; + dnodes[l][n] = 1; + } + } + k += (nnodes[l - 1] + 1) * nnodes[l]; + } + } + + // backwardprop + // output layer dnode initialized to 1. + + for (int n = 0; n < nnodes[nl - 1]; n++) { + if (activation[nl - 1] == 0) { + bnodes[nl - 1][n] = 1; + } else { + bnodes[nl - 1][n] = dnodes[nl - 1][n]; + } + } + + if (nl > 1) { + for (int l = nl - 1; l > 0; l--) { + k -= (nnodes[l - 1] + 1) * nnodes[l]; + for (int n = 0; n < nnodes[l - 1]; n++) { + bnodes[l - 1][n] = 0; + for (int j = 0; j < nnodes[l]; j++) { + bnodes[l - 1][n] += coeffi[k + j * (nnodes[l - 1] + 1) + n + 1] * bnodes[l][j]; + } + if (activation[l - 1] >= 1) { bnodes[l - 1][n] *= dnodes[l - 1][n]; } + } + } + } + + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + data->betas[ii][icoeff] = 0; + for (int j = 0; j < nnodes[0]; j++) { + data->betas[ii][icoeff] += + coeffi[j * ((data->ndescriptors) + 1) + icoeff + 1] * bnodes[0][j]; + } + data->betas[ii][icoeff] = data->betas[ii][icoeff] / scalei[1][icoeff]; + } + + if (data->eflag) { // energy of atom I (E_i) - double etmp = nodes[nl-1][0]; + double etmp = nodes[nl - 1][0]; - data->energy += etmp; - data->eatoms[ii] = etmp; - } - // Deleting the variables + data->energy += etmp; + data->eatoms[ii] = etmp; + } + // Deleting the variables - for (int n=0; nall(FLERR,"compute_gradgrads not implemented"); + error->all(FLERR, "compute_gradgrads not implemented"); } /* ---------------------------------------------------------------------- @@ -391,7 +371,7 @@ void MLIAPModelNN::compute_gradgrads(class MLIAPData * /*data*/) void MLIAPModelNN::compute_force_gradients(class MLIAPData * /*data*/) { - error->all(FLERR,"compute_force_gradients not implemented"); + error->all(FLERR, "compute_force_gradients not implemented"); } /* ---------------------------------------------------------------------- @@ -408,9 +388,9 @@ double MLIAPModelNN::memory_usage() { double bytes = 0; - bytes += (double)nelements*nparams*sizeof(double); // coeffelem - bytes += (double)nelements*2*ndescriptors*sizeof(double); // scale - bytes += (int)nlayers*sizeof(int); // nnodes - bytes += (int)nlayers*sizeof(int); // activation + bytes += (double) nelements * nparams * sizeof(double); // coeffelem + bytes += (double) nelements * 2 * ndescriptors * sizeof(double); // scale + bytes += (int) nlayers * sizeof(int); // nnodes + bytes += (int) nlayers * sizeof(int); // activation return bytes; } diff --git a/src/ML-SNAP/compute_sna_atom.cpp b/src/ML-SNAP/compute_sna_atom.cpp index 3f5c1be3c0..0dc2de4656 100644 --- a/src/ML-SNAP/compute_sna_atom.cpp +++ b/src/ML-SNAP/compute_sna_atom.cpp @@ -32,12 +32,11 @@ using namespace LAMMPS_NS; ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), sna(nullptr), - radelem(nullptr), wjelem(nullptr) + radelem(nullptr), wjelem(nullptr), rinnerelem(nullptr), drinnerelem(nullptr) + { double rmin0, rfac0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; - radelem = nullptr; - wjelem = nullptr; int ntypes = atom->ntypes; int nargmin = 6+2*ntypes; @@ -54,6 +53,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 0; nelements = 1; // offset by 1 to match up with types @@ -133,12 +133,26 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute sna/atom command"); wselfallflag = atoi(arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { + if (iarg+1+2*ntypes > narg) + error->all(FLERR,"Illegal compute sna/atom command"); + switchinnerflag = 1; + iarg++; + memory->create(rinnerelem,ntypes+1,"sna/atom:rinnerelem"); + memory->create(drinnerelem,ntypes+1,"sna/atom:drinnerelem"); + for (int i = 0; i < ntypes; i++) + rinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + iarg += ntypes; + for (int i = 0; i < ntypes; i++) + drinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + iarg += ntypes; } else error->all(FLERR,"Illegal compute sna/atom command"); } snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; size_peratom_cols = ncoeff; @@ -158,6 +172,13 @@ ComputeSNAAtom::~ComputeSNAAtom() memory->destroy(wjelem); memory->destroy(cutsq); delete snaptr; + + if (chemflag) memory->destroy(map); + + if (switchinnerflag) { + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); + } } /* ---------------------------------------------------------------------- */ @@ -268,7 +289,11 @@ void ComputeSNAAtom::compute_peratom() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; - snaptr->element[ninside] = jelem; // element index for multi-element snap + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); + snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); + } + if (chemflag) snaptr->element[ninside] = jelem; ninside++; } } diff --git a/src/ML-SNAP/compute_sna_atom.h b/src/ML-SNAP/compute_sna_atom.h index ca5cd64d47..5e62fe47d9 100644 --- a/src/ML-SNAP/compute_sna_atom.h +++ b/src/ML-SNAP/compute_sna_atom.h @@ -44,6 +44,9 @@ class ComputeSNAAtom : public Compute { double *wjelem; int *map; // map types to [0,nelements) int nelements, chemflag; + int switchinnerflag; + double *rinnerelem; + double *drinnerelem; class SNA *snaptr; double cutmax; int quadraticflag; diff --git a/src/ML-SNAP/compute_snad_atom.cpp b/src/ML-SNAP/compute_snad_atom.cpp index 80076d6afb..25f4b430aa 100644 --- a/src/ML-SNAP/compute_snad_atom.cpp +++ b/src/ML-SNAP/compute_snad_atom.cpp @@ -32,12 +32,10 @@ using namespace LAMMPS_NS; ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snad(nullptr), - radelem(nullptr), wjelem(nullptr) + radelem(nullptr), wjelem(nullptr), rinnerelem(nullptr), drinnerelem(nullptr) { double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; - radelem = nullptr; - wjelem = nullptr; int ntypes = atom->ntypes; int nargmin = 6+2*ntypes; @@ -54,6 +52,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 0; nelements = 1; // process required arguments @@ -131,12 +130,26 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute snad/atom command"); wselfallflag = atoi(arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { + if (iarg+1+2*ntypes > narg) + error->all(FLERR,"Illegal compute snad/atom command"); + switchinnerflag = 1; + iarg++; + memory->create(rinnerelem,ntypes+1,"snad/atom:rinnerelem"); + memory->create(drinnerelem,ntypes+1,"snad/atom:drinnerelem"); + for (int i = 0; i < ntypes; i++) + rinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + iarg += ntypes; + for (int i = 0; i < ntypes; i++) + drinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + iarg += ntypes; } else error->all(FLERR,"Illegal compute snad/atom command"); } snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -160,6 +173,13 @@ ComputeSNADAtom::~ComputeSNADAtom() memory->destroy(wjelem); memory->destroy(cutsq); delete snaptr; + + if (chemflag) memory->destroy(map); + + if (switchinnerflag) { + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); + } } /* ---------------------------------------------------------------------- */ @@ -286,7 +306,11 @@ void ComputeSNADAtom::compute_peratom() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; - snaptr->element[ninside] = jelem; // element index for multi-element snap + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); + snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); + } + if (chemflag) snaptr->element[ninside] = jelem; ninside++; } } @@ -299,8 +323,7 @@ void ComputeSNADAtom::compute_peratom() for (int jj = 0; jj < ninside; jj++) { const int j = snaptr->inside[jj]; - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj], jj, snaptr->element[jj]); + snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate -dBi/dRi, -dBi/dRj diff --git a/src/ML-SNAP/compute_snad_atom.h b/src/ML-SNAP/compute_snad_atom.h index 94dec9b85f..c03458446f 100644 --- a/src/ML-SNAP/compute_snad_atom.h +++ b/src/ML-SNAP/compute_snad_atom.h @@ -46,6 +46,9 @@ class ComputeSNADAtom : public Compute { double *wjelem; int *map; // map types to [0,nelements) int nelements, chemflag; + int switchinnerflag; + double *rinnerelem; + double *drinnerelem; class SNA *snaptr; double cutmax; int quadraticflag; diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp index 4d24616aa5..af8c52c2ee 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -35,7 +35,7 @@ enum{SCALAR,VECTOR,ARRAY}; ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snap(nullptr), snapall(nullptr), snap_peratom(nullptr), radelem(nullptr), wjelem(nullptr), - snaptr(nullptr) + snaptr(nullptr), rinnerelem(nullptr), drinnerelem(nullptr) { array_flag = 1; @@ -43,8 +43,6 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; - radelem = nullptr; - wjelem = nullptr; int ntypes = atom->ntypes; int nargmin = 6+2*ntypes; @@ -61,6 +59,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 0; nelements = 1; // process required arguments @@ -143,12 +142,26 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute snap command"); bikflag = atoi(arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { + if (iarg+1+2*ntypes > narg) + error->all(FLERR,"Illegal compute snap command"); + switchinnerflag = 1; + iarg++; + memory->create(rinnerelem,ntypes+1,"snap:rinnerelem"); + memory->create(drinnerelem,ntypes+1,"snap:drinnerelem"); + for (int i = 0; i < ntypes; i++) + rinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + iarg += ntypes; + for (int i = 0; i < ntypes; i++) + drinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + iarg += ntypes; } else error->all(FLERR,"Illegal compute snap command"); } snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -183,6 +196,11 @@ ComputeSnap::~ComputeSnap() delete snaptr; if (chemflag) memory->destroy(map); + + if (switchinnerflag) { + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); + } } /* ---------------------------------------------------------------------- */ @@ -342,7 +360,11 @@ void ComputeSnap::compute_array() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; - snaptr->element[ninside] = jelem; // element index for multi-element snap + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); + snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); + } + if (chemflag) snaptr->element[ninside] = jelem; ninside++; } } @@ -353,8 +375,7 @@ void ComputeSnap::compute_array() for (int jj = 0; jj < ninside; jj++) { const int j = snaptr->inside[jj]; - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj], jj, snaptr->element[jj]); + snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate dBi/dRi, -dBi/dRj diff --git a/src/ML-SNAP/compute_snap.h b/src/ML-SNAP/compute_snap.h index ccdc0e107e..fa31de747d 100644 --- a/src/ML-SNAP/compute_snap.h +++ b/src/ML-SNAP/compute_snap.h @@ -46,6 +46,9 @@ class ComputeSnap : public Compute { double *wjelem; int *map; // map types to [0,nelements) int nelements, chemflag; + int switchinnerflag; + double *rinnerelem; + double *drinnerelem; class SNA *snaptr; double cutmax; int quadraticflag; diff --git a/src/ML-SNAP/compute_snav_atom.cpp b/src/ML-SNAP/compute_snav_atom.cpp index d0031814ca..9834dc2a0c 100644 --- a/src/ML-SNAP/compute_snav_atom.cpp +++ b/src/ML-SNAP/compute_snav_atom.cpp @@ -31,12 +31,10 @@ using namespace LAMMPS_NS; ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snav(nullptr), - radelem(nullptr), wjelem(nullptr) + radelem(nullptr), wjelem(nullptr), rinnerelem(nullptr), drinnerelem(nullptr) { double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; - radelem = nullptr; - wjelem = nullptr; int ntypes = atom->ntypes; int nargmin = 6+2*ntypes; @@ -53,6 +51,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 0; nelements = 1; // process required arguments @@ -126,12 +125,26 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute snav/atom command"); wselfallflag = atoi(arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { + if (iarg+1+2*ntypes > narg) + error->all(FLERR,"Illegal compute snav/atom command"); + switchinnerflag = 1; + iarg++; + memory->create(rinnerelem,ntypes+1,"snav/atom:rinnerelem"); + memory->create(drinnerelem,ntypes+1,"snav/atom:drinnerelem"); + for (int i = 0; i < ntypes; i++) + rinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + iarg += ntypes; + for (int i = 0; i < ntypes; i++) + drinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + iarg += ntypes; } else error->all(FLERR,"Illegal compute snav/atom command"); } snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -153,8 +166,14 @@ ComputeSNAVAtom::~ComputeSNAVAtom() memory->destroy(radelem); memory->destroy(wjelem); memory->destroy(cutsq); - delete snaptr; + + if (chemflag) memory->destroy(map); + + if (switchinnerflag) { + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); + } } /* ---------------------------------------------------------------------- */ @@ -279,7 +298,11 @@ void ComputeSNAVAtom::compute_peratom() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; - snaptr->element[ninside] = jelem; // element index for multi-element snap + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); + snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); + } + if (chemflag) snaptr->element[ninside] = jelem; ninside++; } } @@ -293,8 +316,7 @@ void ComputeSNAVAtom::compute_peratom() for (int jj = 0; jj < ninside; jj++) { const int j = snaptr->inside[jj]; - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj], jj, snaptr->element[jj]); + snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate -dBi/dRi*Ri, -dBi/dRj*Rj diff --git a/src/ML-SNAP/compute_snav_atom.h b/src/ML-SNAP/compute_snav_atom.h index d61d345408..be442f3a49 100644 --- a/src/ML-SNAP/compute_snav_atom.h +++ b/src/ML-SNAP/compute_snav_atom.h @@ -46,6 +46,9 @@ class ComputeSNAVAtom : public Compute { double *wjelem; int *map; // map types to [0,nelements) int nelements, chemflag; + int switchinnerflag; + double *rinnerelem; + double *drinnerelem; class SNA *snaptr; int quadraticflag; }; diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index eafb27f5ba..e9c6958ceb 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -46,6 +46,8 @@ PairSNAP::PairSNAP(LAMMPS *lmp) : Pair(lmp) radelem = nullptr; wjelem = nullptr; coeffelem = nullptr; + rinnerelem = nullptr; + drinnerelem = nullptr; beta_max = 0; beta = nullptr; @@ -63,6 +65,11 @@ PairSNAP::~PairSNAP() memory->destroy(wjelem); memory->destroy(coeffelem); + if (switchinnerflag) { + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); + } + memory->destroy(beta); memory->destroy(bispectrum); @@ -151,7 +158,11 @@ void PairSNAP::compute(int eflag, int vflag) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; - snaptr->element[ninside] = jelem; + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); + snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + } + if (chemflag) snaptr->element[ninside] = jelem; ninside++; } } @@ -172,12 +183,7 @@ void PairSNAP::compute(int eflag, int vflag) for (int jj = 0; jj < ninside; jj++) { int j = snaptr->inside[jj]; - if (chemflag) - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, snaptr->element[jj]); - else - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, 0); + snaptr->compute_duidrj(jj); snaptr->compute_deidrj(fij); @@ -325,7 +331,11 @@ void PairSNAP::compute_bispectrum() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; - snaptr->element[ninside] = jelem; + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); + snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + } + if (chemflag) snaptr->element[ninside] = jelem; ninside++; } } @@ -403,7 +413,8 @@ void PairSNAP::coeff(int narg, char **arg) snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, + nelements, switchinnerflag); if (ncoeff != snaptr->ncoeff) { if (comm->me == 0) @@ -508,9 +519,13 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) memory->destroy(radelem); memory->destroy(wjelem); memory->destroy(coeffelem); + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); memory->create(radelem,nelements,"pair:radelem"); memory->create(wjelem,nelements,"pair:wjelem"); memory->create(coeffelem,nelements,ncoeffall,"pair:coeffelem"); + memory->create(rinnerelem,nelements,"pair:rinnerelem"); + memory->create(drinnerelem,nelements,"pair:drinnerelem"); // initialize checklist for all required nelements @@ -626,9 +641,15 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 0; chunksize = 32768; parallel_thresh = 8192; + // set local input checks + + int rinnerflag = 0; + int drinnerflag = 0; + // open SNAP parameter file on proc 0 FILE *fpparam; @@ -652,7 +673,8 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) if (eof) break; MPI_Bcast(line,MAXLINE,MPI_CHAR,0,world); - // strip comment, skip line if blank + // words = ptrs to all words in line + // strip single and double quotes from words std::vector words; try { @@ -662,43 +684,83 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) } if (words.size() == 0) continue; - if (words.size() != 2) + + if (words.size() < 2) error->all(FLERR,"Incorrect format in SNAP parameter file"); auto keywd = words[0]; auto keyval = words[1]; - if (comm->me == 0) - utils::logmesg(lmp,"SNAP keyword {} {}\n",keywd,keyval); + // check for keywords with one value per element - if (keywd == "rcutfac") { - rcutfac = utils::numeric(FLERR,keyval,false,lmp); - rcutfacflag = 1; - } else if (keywd == "twojmax") { - twojmax = utils::inumeric(FLERR,keyval,false,lmp); - twojmaxflag = 1; - } else if (keywd == "rfac0") - rfac0 = utils::numeric(FLERR,keyval,false,lmp); - else if (keywd == "rmin0") - rmin0 = utils::numeric(FLERR,keyval,false,lmp); - else if (keywd == "switchflag") - switchflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (keywd == "bzeroflag") - bzeroflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (keywd == "quadraticflag") - quadraticflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (keywd == "chemflag") - chemflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (keywd == "bnormflag") - bnormflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (keywd == "wselfallflag") - wselfallflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (keywd == "chunksize") - chunksize = utils::inumeric(FLERR,keyval,false,lmp); - else if (keywd == "parallelthresh") - parallel_thresh = utils::inumeric(FLERR,keyval,false,lmp); - else - error->all(FLERR,"Unknown parameter '{}' in SNAP parameter file", keywd); + if (keywd == "rinner" || keywd == "drinner") { + + if (nwords != nelements+1) + error->all(FLERR,"Incorrect SNAP parameter file"); + + if (comm->me == 0) + utils::logmesg(lmp,"SNAP keyword {} {} ... \n", keywd, keyval); + + int iword = 1; + + if (keywd == "rinner") { + keyval = words[iword]; + for (int ielem = 0; ielem < nelements; ielem++) { + printf("rinnerelem = %p ielem = %d nelements = %d iword = %d nwords = %d\n",rinnerelem, ielem, nelements, iword, nwords); + rinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + iword++; + } + rinnerflag = 1; + } else if (keywd == "drinner") { + keyval = words[iword]; + for (int ielem = 0; ielem < nelements; ielem++) { + drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + iword++; + } + drinnerflag = 1; + } + + } else { + + // all other keywords take one value + + if (nwords != 2) + error->all(FLERR,"Incorrect SNAP parameter file"); + + if (comm->me == 0) + utils::logmesg(lmp,"SNAP keyword {} {}\n",keywd,keyval); + + if (keywd == "rcutfac") { + rcutfac = utils::numeric(FLERR,keyval,false,lmp); + rcutfacflag = 1; + } else if (keywd == "twojmax") { + twojmax = utils::inumeric(FLERR,keyval,false,lmp); + twojmaxflag = 1; + } else if (keywd == "rfac0") + rfac0 = utils::numeric(FLERR,keyval,false,lmp); + else if (keywd == "rmin0") + rmin0 = utils::numeric(FLERR,keyval,false,lmp); + else if (keywd == "switchflag") + switchflag = utils::inumeric(FLERR,keyval,false,lmp); + else if (keywd == "bzeroflag") + bzeroflag = utils::inumeric(FLERR,keyval,false,lmp); + else if (keywd == "quadraticflag") + quadraticflag = utils::inumeric(FLERR,keyval,false,lmp); + else if (keywd == "chemflag") + chemflag = utils::inumeric(FLERR,keyval,false,lmp); + else if (keywd == "bnormflag") + bnormflag = utils::inumeric(FLERR,keyval,false,lmp); + else if (keywd == "wselfallflag") + wselfallflag = utils::inumeric(FLERR,keyval,false,lmp); + else if (keywd == "switchinnerflag") + switchinnerflag = utils::inumeric(FLERR,keyval,false,lmp); + else if (keywd == "chunksize") + chunksize = utils::inumeric(FLERR,keyval,false,lmp); + else if (keywd == "parallelthresh") + parallel_thresh = utils::inumeric(FLERR,keyval,false,lmp); + else + error->all(FLERR,"Unknown parameter '{}' in SNAP parameter file", keywd); + } } if (rcutfacflag == 0 || twojmaxflag == 0) @@ -707,6 +769,11 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) if (chemflag && nelemtmp != nelements) error->all(FLERR,"Incorrect SNAP parameter file"); + if (switchinnerflag && !(rinnerflag && drinnerflag)) + error->all(FLERR,"Incorrect SNAP parameter file"); + + if (!switchinnerflag && (rinnerflag || drinnerflag)) + error->all(FLERR,"Incorrect SNAP parameter file"); } /* ---------------------------------------------------------------------- diff --git a/src/ML-SNAP/pair_snap.h b/src/ML-SNAP/pair_snap.h index 26640aed1a..77ceeabd4f 100644 --- a/src/ML-SNAP/pair_snap.h +++ b/src/ML-SNAP/pair_snap.h @@ -59,6 +59,9 @@ class PairSNAP : public Pair { double **scale; // for thermodynamic integration int twojmax, switchflag, bzeroflag, bnormflag; int chemflag, wselfallflag; + int switchinnerflag; // inner cutoff switch + double *rinnerelem; // element inner cutoff + double *drinnerelem; // element inner cutoff range int chunksize,parallel_thresh; double rfac0, rmin0, wj1, wj2; int rcutfacflag, twojmaxflag; // flags for required parameters diff --git a/src/ML-SNAP/sna.cpp b/src/ML-SNAP/sna.cpp index 8a7ffa9f45..b33ad8a491 100644 --- a/src/ML-SNAP/sna.cpp +++ b/src/ML-SNAP/sna.cpp @@ -103,7 +103,7 @@ using namespace MathSpecial; j = |j1-j2|, |j1-j2|+2,...,j1+j2-2,j1+j2 [1] Albert Bartok-Partay, "Gaussian Approximation..." - Doctoral Thesis, Cambrindge University, (2009) + Doctoral Thesis, Cambridge University, (2009) [2] D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii, "Quantum Theory of Angular Momentum," World Scientific (1988) @@ -112,13 +112,15 @@ using namespace MathSpecial; SNA::SNA(LAMMPS* lmp, double rfac0_in, int twojmax_in, double rmin0_in, int switch_flag_in, int bzero_flag_in, - int chem_flag_in, int bnorm_flag_in, int wselfall_flag_in, int nelements_in) : Pointers(lmp) + int chem_flag_in, int bnorm_flag_in, int wselfall_flag_in, + int nelements_in, int switch_inner_flag_in) : Pointers(lmp) { wself = 1.0; rfac0 = rfac0_in; rmin0 = rmin0_in; switch_flag = switch_flag_in; + switch_inner_flag = switch_inner_flag_in; bzero_flag = bzero_flag_in; chem_flag = chem_flag_in; bnorm_flag = bnorm_flag_in; @@ -141,6 +143,8 @@ SNA::SNA(LAMMPS* lmp, double rfac0_in, int twojmax_in, inside = nullptr; wj = nullptr; rcutij = nullptr; + rinnerij = nullptr; + drinnerij = nullptr; element = nullptr; nmax = 0; idxz = nullptr; @@ -170,7 +174,11 @@ SNA::~SNA() memory->destroy(inside); memory->destroy(wj); memory->destroy(rcutij); - memory->destroy(element); + if (switch_inner_flag) { + memory->destroy(rinnerij); + memory->destroy(drinnerij); + } + if (chem_flag) memory->destroy(element); memory->destroy(ulist_r_ij); memory->destroy(ulist_i_ij); delete[] idxz; @@ -307,7 +315,6 @@ void SNA::init() init_rootpqarray(); } - void SNA::grow_rij(int newnmax) { if (newnmax <= nmax) return; @@ -318,14 +325,22 @@ void SNA::grow_rij(int newnmax) memory->destroy(inside); memory->destroy(wj); memory->destroy(rcutij); - memory->destroy(element); + if (switch_inner_flag) { + memory->destroy(rinnerij); + memory->destroy(drinnerij); + } + if (chem_flag) memory->destroy(element); memory->destroy(ulist_r_ij); memory->destroy(ulist_i_ij); memory->create(rij, nmax, 3, "pair:rij"); memory->create(inside, nmax, "pair:inside"); memory->create(wj, nmax, "pair:wj"); memory->create(rcutij, nmax, "pair:rcutij"); - memory->create(element, nmax, "sna:element"); + if (switch_inner_flag) { + memory->create(rinnerij, nmax, "pair:rinnerij"); + memory->create(drinnerij, nmax, "pair:drinnerij"); + } + if (chem_flag) memory->create(element, nmax, "sna:element"); memory->create(ulist_r_ij, nmax, idxu_max, "sna:ulist_ij"); memory->create(ulist_i_ij, nmax, idxu_max, "sna:ulist_ij"); } @@ -358,10 +373,7 @@ void SNA::compute_ui(int jnum, int ielem) z0 = r / tan(theta0); compute_uarray(x, y, z, z0, r, j); - if (chem_flag) - add_uarraytot(r, wj[j], rcutij[j], j, element[j]); - else - add_uarraytot(r, wj[j], rcutij[j], j, 0); + add_uarraytot(r, j); } } @@ -951,14 +963,15 @@ void SNA::compute_dbidrj() calculate derivative of Ui w.r.t. atom j ------------------------------------------------------------------------- */ -void SNA::compute_duidrj(double* rij, double wj, double rcut, int jj, int jelem) +void SNA::compute_duidrj(int jj) { double rsq, r, x, y, z, z0, theta0, cs, sn; double dz0dr; + double rcut = rcutij[jj]; - x = rij[0]; - y = rij[1]; - z = rij[2]; + x = rij[jj][0]; + y = rij[jj][1]; + z = rij[jj][2]; rsq = x * x + y * y + z * z; r = sqrt(rsq); double rscale0 = rfac0 * MY_PI / (rcut - rmin0); @@ -968,8 +981,10 @@ void SNA::compute_duidrj(double* rij, double wj, double rcut, int jj, int jelem) z0 = r * cs / sn; dz0dr = z0 / r - (r*rscale0) * (rsq + z0 * z0) / rsq; - elem_duarray = jelem; - compute_duarray(x, y, z, z0, r, dz0dr, wj, rcut, jj); + if (chem_flag) elem_duarray = element[jj]; + else elem_duarray = 0; + + compute_duarray(x, y, z, z0, r, dz0dr, wj[jj], rcut, jj); } /* ---------------------------------------------------------------------- */ @@ -998,13 +1013,19 @@ void SNA::zero_uarraytot(int ielem) add Wigner U-functions for one neighbor to the total ------------------------------------------------------------------------- */ -void SNA::add_uarraytot(double r, double wj, double rcut, int jj, int jelem) +void SNA::add_uarraytot(double r, int jj) { double sfac; + int jelem; - sfac = compute_sfac(r, rcut); + sfac = compute_sfac(r, rcutij[jj]); + sfac *= wj[jj]; - sfac *= wj; + if (switch_inner_flag) + sfac *= compute_sfac_inner(r, rinnerij[jj], drinnerij[jj]); + + if (chem_flag) jelem = element[jj]; + else jelem = 0; double* ulist_r = ulist_r_ij[jj]; double* ulist_i = ulist_i_ij[jj]; @@ -1248,6 +1269,12 @@ void SNA::compute_duarray(double x, double y, double z, double sfac = compute_sfac(r, rcut); double dsfac = compute_dsfac(r, rcut); + if (switch_inner_flag) { + double sfac_inner = compute_sfac_inner(r, rinnerij[jj], drinnerij[jj]); + dsfac = dsfac*sfac_inner + sfac*compute_dsfac_inner(r, rinnerij[jj], drinnerij[jj]); + sfac *= sfac_inner; + } + sfac *= wj; dsfac *= wj; for (int j = 0; j <= twojmax; j++) { @@ -1310,6 +1337,11 @@ double SNA::memory_usage() bytes += (double)nmax * sizeof(int); // inside bytes += (double)nmax * sizeof(double); // wj bytes += (double)nmax * sizeof(double); // rcutij + if (switch_inner_flag) { + bytes += (double)nmax * sizeof(double); // rinnerij + bytes += (double)nmax * sizeof(double); // drinnerij + } + if (chem_flag) bytes += (double)nmax * sizeof(int); // element return bytes; } @@ -1515,31 +1547,59 @@ void SNA::compute_ncoeff() double SNA::compute_sfac(double r, double rcut) { - if (switch_flag == 0) return 1.0; - if (switch_flag == 1) { - if (r <= rmin0) return 1.0; - else if (r > rcut) return 0.0; - else { - double rcutfac = MY_PI / (rcut - rmin0); - return 0.5 * (cos((r - rmin0) * rcutfac) + 1.0); - } + double sfac; + if (switch_flag == 0) sfac = 1.0; + else if (r <= rmin0) sfac = 1.0; + else if (r > rcut) sfac = 0.0; + else { + double rcutfac = MY_PI / (rcut - rmin0); + sfac = 0.5 * (cos((r - rmin0) * rcutfac) + 1.0); } - return 0.0; + return sfac; } /* ---------------------------------------------------------------------- */ double SNA::compute_dsfac(double r, double rcut) { - if (switch_flag == 0) return 0.0; - if (switch_flag == 1) { - if (r <= rmin0) return 0.0; - else if (r > rcut) return 0.0; - else { - double rcutfac = MY_PI / (rcut - rmin0); - return -0.5 * sin((r - rmin0) * rcutfac) * rcutfac; - } + double dsfac; + if (switch_flag == 0) dsfac = 0.0; + else if (r <= rmin0) dsfac = 0.0; + else if (r > rcut) dsfac = 0.0; + else { + double rcutfac = MY_PI / (rcut - rmin0); + dsfac = -0.5 * sin((r - rmin0) * rcutfac) * rcutfac; } - return 0.0; + return dsfac; +} + +/* ---------------------------------------------------------------------- */ + +double SNA::compute_sfac_inner(double r, double rinner, double drinner) +{ + double sfac; + if (switch_inner_flag == 0) sfac = 1.0; + else if (r >= rinner + drinner) sfac = 1.0; + else if (r <= rinner) sfac = 0.0; + else { + double rcutfac = MY_PI / drinner; + sfac = 0.5 * (1.0 - cos((r - rinner) * rcutfac)); + } + return sfac; +} + +/* ---------------------------------------------------------------------- */ + +double SNA::compute_dsfac_inner(double r, double rinner, double drinner) +{ + double dsfac; + if (switch_inner_flag == 0) dsfac = 0.0; + else if (r >= rinner + drinner) dsfac = 0.0; + else if (r <= rinner) dsfac = 0.0; + else { + double rcutfac = MY_PI / drinner; + dsfac = 0.5 * sin((r - rinner) * rcutfac) * rcutfac; + } + return dsfac; } diff --git a/src/ML-SNAP/sna.h b/src/ML-SNAP/sna.h index 3b949d564c..03d77d52b4 100644 --- a/src/ML-SNAP/sna.h +++ b/src/ML-SNAP/sna.h @@ -33,7 +33,7 @@ struct SNA_BINDICES { class SNA : protected Pointers { public: - SNA(LAMMPS *, double, int, double, int, int, int, int, int, int); + SNA(LAMMPS *, double, int, double, int, int, int, int, int, int, int); SNA(LAMMPS *lmp) : Pointers(lmp){}; ~SNA() override; @@ -53,26 +53,39 @@ class SNA : protected Pointers { // functions for derivatives - void compute_duidrj(double *, double, double, int, int); + void compute_duidrj(int); void compute_dbidrj(); void compute_deidrj(double *); double compute_sfac(double, double); double compute_dsfac(double, double); - double *blist; - double **dblist; - double **rij; - int *inside; - double *wj; - double *rcutij; - int *element; // index on [0,nelements) - int nmax; + double compute_sfac_inner(double, double, double); + double compute_dsfac_inner(double, double, double); - void grow_rij(int); + // public bispectrum data int twojmax; - double *ylist_r, *ylist_i; - int idxcg_max, idxu_max, idxz_max, idxb_max; + double *blist; + double **dblist; + + // short neighbor list data + + void grow_rij(int); + int nmax; // allocated size of short lists + + double **rij; // short rij list + int *inside; // short neighbor list + double *wj; // short weight list + double *rcutij; // short cutoff list + + // only allocated for switch_inner_flag=1 + + double *rinnerij; // short inner cutoff list + double *drinnerij;// short inner range list + + // only allocated for chem_flag=1 + + int *element; // short element list [0,nelements) private: double rmin0, rfac0; @@ -87,7 +100,7 @@ class SNA : protected Pointers { int ***idxcg_block; double *ulisttot_r, *ulisttot_i; - double **ulist_r_ij, **ulist_i_ij; + double **ulist_r_ij, **ulist_i_ij; // short u list int *idxu_block; double *zlist_r, *zlist_i; @@ -98,23 +111,32 @@ class SNA : protected Pointers { double **dulist_r, **dulist_i; int elem_duarray; // element of j in derivative + double *ylist_r, *ylist_i; + int idxcg_max, idxu_max, idxz_max, idxb_max; + void create_twojmax_arrays(); void destroy_twojmax_arrays(); void init_clebsch_gordan(); void print_clebsch_gordan(); void init_rootpqarray(); void zero_uarraytot(int); - void add_uarraytot(double, double, double, int, int); + void add_uarraytot(double, int); void compute_uarray(double, double, double, double, double, int); double deltacg(int, int, int); void compute_ncoeff(); - void compute_duarray(double, double, double, double, double, double, double, double, int); + void compute_duarray(double, double, double, double, double, double, + double, double, int); // Sets the style for the switching function // 0 = none // 1 = cosine int switch_flag; + // Sets the style for the inner switching function + // 0 = none + // 1 = cosine + int switch_inner_flag; + // Self-weight double wself;