ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3173 from athomps/snap-inner

Browse Source 
Inner cutoff switching function for SNAP descriptors and potentials
This commit is contained in:
Axel Kohlmeyer
2022-03-16 12:43:23 -04:00
committed by GitHub
parent 4894263215 8e585620c3
commit ba19090631
24 changed files with 753 additions and 461 deletions
Download Patch File Download Diff File Expand all files Collapse all files

26
doc/src/compute_sna_atom.rst
View File

@ -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

14
doc/src/pair_snap.rst
View File

@ -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

2
doc/utils/sphinx-config/false_positives.txt
View File

@ -1928,6 +1928,7 @@ mbt
MBytes
mc
McLachlan
mcmoves
md
mdf
mdi
@ -2228,6 +2229,7 @@ Nelement
Nelements
nelems
nemd
netapp
netcdf
netstat
Nettleton

2
potentials/InP_JCPA2020.mliap
View File

@ -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.

2
potentials/InP_JCPA2020.mliap.descriptor
View File

@ -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

2
potentials/InP_JCPA2020.mliap.model
View File

@ -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

2
potentials/InP_JCPA2020.snap
View File

@ -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.

2
potentials/InP_JCPA2020.snapcoeff
View File

@ -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

2
potentials/InP_JCPA2020.snapparam
View File

@ -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

263
src/ML-IAP/mliap_descriptor_snap.cpp
View File

@ -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 <cmath>
#include <cstring>
@ -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;
}

4
src/ML-IAP/mliap_descriptor_snap.h
View File

@ -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

396
src/ML-IAP/mliap_model_nn.cpp
View File

@ -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; l<nl; ++l) {
nodes[l] = new double[nnodes[l]];
dnodes[l] = new double[nnodes[l]];
bnodes[l] = new double[nnodes[l]];
}
nodes = new double *[nl];
dnodes = new double *[nl];
bnodes = new double *[nl];
for (int l = 0; l < nl; ++l) {
nodes[l] = new double[nnodes[l]];
dnodes[l] = new double[nnodes[l]];
bnodes[l] = new double[nnodes[l]];
}
// forwardprop
// input - hidden1
for (int n=0; n < nnodes[0]; n++) {
nodes[0][n] = 0;
for (int icoeff = 0; icoeff < data->ndescriptors; 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<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];
}
}
}
}
// 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; 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];
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; n<nl; n++) {
delete[] nodes[n];
delete[] dnodes[n];
delete[] bnodes[n];
}
delete[] nodes;
delete[] dnodes;
delete[] bnodes;
for (int n = 0; n < nl; n++) {
delete[] nodes[n];
delete[] dnodes[n];
delete[] bnodes[n];
}
delete[] nodes;
delete[] dnodes;
delete[] bnodes;
}
}
@ -381,7 +361,7 @@ void MLIAPModelNN::compute_gradients(MLIAPData* data)
void MLIAPModelNN::compute_gradgrads(class MLIAPData * /*data*/)
{
error->all(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;
}

35
src/ML-SNAP/compute_sna_atom.cpp
View File

@ -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++;
}
}

3
src/ML-SNAP/compute_sna_atom.h
View File

@ -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;

37
src/ML-SNAP/compute_snad_atom.cpp
View File

@ -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

3
src/ML-SNAP/compute_snad_atom.h
View File

@ -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;

35
src/ML-SNAP/compute_snap.cpp
View File

@ -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

3
src/ML-SNAP/compute_snap.h
View File

@ -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;

38
src/ML-SNAP/compute_snav_atom.cpp
View File

@ -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

3
src/ML-SNAP/compute_snav_atom.h
View File

@ -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;
};

149
src/ML-SNAP/pair_snap.cpp
View File

@ -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<std::string> 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");
}
/* ----------------------------------------------------------------------

3
src/ML-SNAP/pair_snap.h
View File

@ -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

134
src/ML-SNAP/sna.cpp
View File

@ -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;
}

54
src/ML-SNAP/sna.h
View File

@ -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;