From adef6bd10cfe2e325a33e22bd8f45d90ba6a8f4d Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Thu, 10 Mar 2022 13:58:14 -0700 Subject: [PATCH 01/16] Added first draft of inner switching functions --- src/ML-SNAP/sna.cpp | 79 ++++++++++++++++++++++++++++++++++----------- src/ML-SNAP/sna.h | 8 +++++ 2 files changed, 68 insertions(+), 19 deletions(-) diff --git a/src/ML-SNAP/sna.cpp b/src/ML-SNAP/sna.cpp index 8a7ffa9f45..9a4c56cc8c 100644 --- a/src/ML-SNAP/sna.cpp +++ b/src/ML-SNAP/sna.cpp @@ -119,6 +119,7 @@ SNA::SNA(LAMMPS* lmp, double rfac0_in, int twojmax_in, rfac0 = rfac0_in; rmin0 = rmin0_in; switch_flag = switch_flag_in; + switch_inner_flag = 0; bzero_flag = bzero_flag_in; chem_flag = chem_flag_in; bnorm_flag = bnorm_flag_in; @@ -1001,11 +1002,15 @@ void SNA::zero_uarraytot(int ielem) void SNA::add_uarraytot(double r, double wj, double rcut, int jj, int jelem) { double sfac; + double rinner = 0.0; + double drinner = 1.0; sfac = compute_sfac(r, rcut); - sfac *= wj; + if (switch_inner_flag) + sfac *= compute_sfac_inner(r, rinner, drinner); + double* ulist_r = ulist_r_ij[jj]; double* ulist_i = ulist_i_ij[jj]; @@ -1248,6 +1253,14 @@ 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 rinner = 0.0; + double drinner = 1.0; + + sfac *= compute_sfac_inner(r, rinner, drinner); + dsfac *= compute_dsfac_inner(r, rinner, drinner); + } + sfac *= wj; dsfac *= wj; for (int j = 0; j <= twojmax; j++) { @@ -1515,31 +1528,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..d5bbd3f304 100644 --- a/src/ML-SNAP/sna.h +++ b/src/ML-SNAP/sna.h @@ -59,6 +59,9 @@ class SNA : protected Pointers { double compute_sfac(double, double); double compute_dsfac(double, double); + double compute_sfac_inner(double, double, double); + double compute_dsfac_inner(double, double, double); + double *blist; double **dblist; double **rij; @@ -115,6 +118,11 @@ class SNA : protected Pointers { // 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; From 11d3b13e31f48533c110ee167cf265d42c380d33 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Thu, 10 Mar 2022 18:40:59 -0700 Subject: [PATCH 02/16] Added a lot of boilerplate for inner data, not working yet --- src/ML-IAP/mliap_descriptor_snap.cpp | 22 ++++++---- src/ML-IAP/mliap_descriptor_snap.h | 1 + src/ML-SNAP/compute_sna_atom.cpp | 4 +- src/ML-SNAP/compute_sna_atom.h | 1 + src/ML-SNAP/compute_snad_atom.cpp | 7 +++- src/ML-SNAP/compute_snad_atom.h | 1 + src/ML-SNAP/compute_snap.cpp | 7 +++- src/ML-SNAP/compute_snap.h | 1 + src/ML-SNAP/compute_snav_atom.cpp | 7 +++- src/ML-SNAP/compute_snav_atom.h | 1 + src/ML-SNAP/pair_snap.cpp | 11 +++-- src/ML-SNAP/pair_snap.h | 1 + src/ML-SNAP/sna.cpp | 62 ++++++++++++++++++---------- src/ML-SNAP/sna.h | 25 ++++++----- 14 files changed, 101 insertions(+), 50 deletions(-) diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index 5475117a99..6571637681 100644 --- a/src/ML-IAP/mliap_descriptor_snap.cpp +++ b/src/ML-IAP/mliap_descriptor_snap.cpp @@ -46,7 +46,8 @@ MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename): snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, + nelements, switchinnerflag); ndescriptors = snaptr->ncoeff; } @@ -162,10 +163,12 @@ void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data) int j = snaptr->inside[jj]; if (chemflag) snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, snaptr->element[jj]); + snaptr->rcutij[jj], jj, snaptr->element[jj], + snaptr->rinnerij[jj], snaptr->drinnerij[jj]); else snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, 0); + snaptr->rcutij[jj], jj, 0, + snaptr->rinnerij[jj], snaptr->drinnerij[jj]); snaptr->compute_deidrj(fij); @@ -236,10 +239,12 @@ void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data) if (chemflag) snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, snaptr->element[jj]); + snaptr->rcutij[jj], jj, snaptr->element[jj], + snaptr->rinnerij[jj], snaptr->drinnerij[jj]); else snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, 0); + snaptr->rcutij[jj], jj, 0, + snaptr->rinnerij[jj], snaptr->drinnerij[jj]); snaptr->compute_dbidrj(); @@ -309,10 +314,12 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data) 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]); + snaptr->rcutij[jj], jj, snaptr->element[jj], + snaptr->rinnerij[jj], snaptr->drinnerij[jj]); else snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, 0); + snaptr->rcutij[jj], jj, 0, + snaptr->rinnerij[jj], snaptr->drinnerij[jj]); snaptr->compute_dbidrj(); @@ -361,6 +368,7 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 1; for (int i = 0; i < nelements; i++) delete[] elements[i]; delete[] elements; diff --git a/src/ML-IAP/mliap_descriptor_snap.h b/src/ML-IAP/mliap_descriptor_snap.h index 9098015a4f..8c8a8f16f4 100644 --- a/src/ML-IAP/mliap_descriptor_snap.h +++ b/src/ML-IAP/mliap_descriptor_snap.h @@ -39,6 +39,7 @@ class MLIAPDescriptorSNAP : public MLIAPDescriptor { int twojmax, switchflag, bzeroflag; int chemflag, bnormflag, wselfallflag; + int switchinnerflag; double rfac0, rmin0; }; diff --git a/src/ML-SNAP/compute_sna_atom.cpp b/src/ML-SNAP/compute_sna_atom.cpp index 3f5c1be3c0..5dd337286b 100644 --- a/src/ML-SNAP/compute_sna_atom.cpp +++ b/src/ML-SNAP/compute_sna_atom.cpp @@ -54,6 +54,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 1; nelements = 1; // offset by 1 to match up with types @@ -138,7 +139,8 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : 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; diff --git a/src/ML-SNAP/compute_sna_atom.h b/src/ML-SNAP/compute_sna_atom.h index ca5cd64d47..e741e02bc2 100644 --- a/src/ML-SNAP/compute_sna_atom.h +++ b/src/ML-SNAP/compute_sna_atom.h @@ -44,6 +44,7 @@ class ComputeSNAAtom : public Compute { double *wjelem; int *map; // map types to [0,nelements) int nelements, chemflag; + int switchinnerflag; 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..a94555faa2 100644 --- a/src/ML-SNAP/compute_snad_atom.cpp +++ b/src/ML-SNAP/compute_snad_atom.cpp @@ -54,6 +54,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 1; nelements = 1; // process required arguments @@ -136,7 +137,8 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -300,7 +302,8 @@ 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->rcutij[jj], jj, snaptr->element[jj], + snaptr->rinnerij[jj], snaptr->drinnerij[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..7727a4653b 100644 --- a/src/ML-SNAP/compute_snad_atom.h +++ b/src/ML-SNAP/compute_snad_atom.h @@ -46,6 +46,7 @@ class ComputeSNADAtom : public Compute { double *wjelem; int *map; // map types to [0,nelements) int nelements, chemflag; + int switchinnerflag; 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..6cf7a38235 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -61,6 +61,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 1; nelements = 1; // process required arguments @@ -148,7 +149,8 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -354,7 +356,8 @@ 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->rcutij[jj], jj, snaptr->element[jj], + snaptr->rinnerij[jj], snaptr->drinnerij[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..9d928a8166 100644 --- a/src/ML-SNAP/compute_snap.h +++ b/src/ML-SNAP/compute_snap.h @@ -46,6 +46,7 @@ class ComputeSnap : public Compute { double *wjelem; int *map; // map types to [0,nelements) int nelements, chemflag; + int switchinnerflag; 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..dd49a90b90 100644 --- a/src/ML-SNAP/compute_snav_atom.cpp +++ b/src/ML-SNAP/compute_snav_atom.cpp @@ -53,6 +53,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 1; nelements = 1; // process required arguments @@ -131,7 +132,8 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -294,7 +296,8 @@ void ComputeSNAVAtom::compute_peratom() const int j = snaptr->inside[jj]; snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj], jj, snaptr->element[jj]); + snaptr->rcutij[jj], jj, snaptr->element[jj], + snaptr->rinnerij[jj], snaptr->drinnerij[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..ea5ea50086 100644 --- a/src/ML-SNAP/compute_snav_atom.h +++ b/src/ML-SNAP/compute_snav_atom.h @@ -46,6 +46,7 @@ class ComputeSNAVAtom : public Compute { double *wjelem; int *map; // map types to [0,nelements) int nelements, chemflag; + int switchinnerflag; class SNA *snaptr; int quadraticflag; }; diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index eafb27f5ba..e263197eb2 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -174,10 +174,13 @@ void PairSNAP::compute(int eflag, int vflag) int j = snaptr->inside[jj]; if (chemflag) snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, snaptr->element[jj]); + snaptr->rcutij[jj], jj, snaptr->element[jj], + snaptr->rinnerij[jj], snaptr->drinnerij[jj]); else snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj],jj, 0); + snaptr->rcutij[jj], jj, 0, + // snaptr->rinnerij[jj], snaptr->drinnerij[jj],jj); + 0.0, 1.0); snaptr->compute_deidrj(fij); @@ -403,7 +406,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) @@ -626,6 +630,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) chemflag = 0; bnormflag = 0; wselfallflag = 0; + switchinnerflag = 1; chunksize = 32768; parallel_thresh = 8192; diff --git a/src/ML-SNAP/pair_snap.h b/src/ML-SNAP/pair_snap.h index 26640aed1a..5cab36d1c4 100644 --- a/src/ML-SNAP/pair_snap.h +++ b/src/ML-SNAP/pair_snap.h @@ -59,6 +59,7 @@ class PairSNAP : public Pair { double **scale; // for thermodynamic integration int twojmax, switchflag, bzeroflag, bnormflag; int chemflag, wselfallflag; + int switchinnerflag; 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 9a4c56cc8c..58a1597a4a 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,14 +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 = 0; + switch_inner_flag = switch_inner_flag_in; bzero_flag = bzero_flag_in; chem_flag = chem_flag_in; bnorm_flag = bnorm_flag_in; @@ -142,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; @@ -171,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; @@ -308,7 +315,6 @@ void SNA::init() init_rootpqarray(); } - void SNA::grow_rij(int newnmax) { if (newnmax <= nmax) return; @@ -319,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"); } @@ -360,9 +374,9 @@ void SNA::compute_ui(int jnum, int ielem) compute_uarray(x, y, z, z0, r, j); if (chem_flag) - add_uarraytot(r, wj[j], rcutij[j], j, element[j]); + add_uarraytot(r, wj[j], rcutij[j], j, element[j], rinnerij[j], drinnerij[j]); else - add_uarraytot(r, wj[j], rcutij[j], j, 0); + add_uarraytot(r, wj[j], rcutij[j], j, 0, rinnerij[j], drinnerij[j]); } } @@ -952,7 +966,8 @@ 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(double* rij, double wj, double rcut, int jj, + int jelem, double rinner, double drinner) { double rsq, r, x, y, z, z0, theta0, cs, sn; double dz0dr; @@ -970,7 +985,7 @@ void SNA::compute_duidrj(double* rij, double wj, double rcut, int jj, int jelem) dz0dr = z0 / r - (r*rscale0) * (rsq + z0 * z0) / rsq; elem_duarray = jelem; - compute_duarray(x, y, z, z0, r, dz0dr, wj, rcut, jj); + compute_duarray(x, y, z, z0, r, dz0dr, wj, rcut, jj, rinner, drinner); } /* ---------------------------------------------------------------------- */ @@ -999,18 +1014,17 @@ 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, double wj, double rcut, int jj, int jelem, + double rinner, double drinner) { double sfac; - double rinner = 0.0; - double drinner = 1.0; sfac = compute_sfac(r, rcut); sfac *= wj; if (switch_inner_flag) sfac *= compute_sfac_inner(r, rinner, drinner); - + double* ulist_r = ulist_r_ij[jj]; double* ulist_i = ulist_i_ij[jj]; @@ -1123,7 +1137,8 @@ void SNA::compute_uarray(double x, double y, double z, void SNA::compute_duarray(double x, double y, double z, double z0, double r, double dz0dr, - double wj, double rcut, int jj) + double wj, double rcut, int jj, + double rinner, double drinner) { double r0inv; double a_r, a_i, b_r, b_i; @@ -1254,11 +1269,9 @@ void SNA::compute_duarray(double x, double y, double z, double dsfac = compute_dsfac(r, rcut); if (switch_inner_flag) { - double rinner = 0.0; - double drinner = 1.0; - - sfac *= compute_sfac_inner(r, rinner, drinner); - dsfac *= compute_dsfac_inner(r, rinner, drinner); + double sfac_inner = compute_sfac_inner(r, rinner, drinner); + dsfac = dsfac*sfac_inner + sfac*compute_dsfac_inner(r, rinner, drinner); + sfac *= sfac_inner; } sfac *= wj; @@ -1323,7 +1336,12 @@ 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; } /* ---------------------------------------------------------------------- */ diff --git a/src/ML-SNAP/sna.h b/src/ML-SNAP/sna.h index d5bbd3f304..8c631abccd 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,7 +53,7 @@ class SNA : protected Pointers { // functions for derivatives - void compute_duidrj(double *, double, double, int, int); + void compute_duidrj(double *, double, double, int, int, double, double); void compute_dbidrj(); void compute_deidrj(double *); double compute_sfac(double, double); @@ -64,12 +64,14 @@ class SNA : protected Pointers { double *blist; double **dblist; - double **rij; - int *inside; - double *wj; - double *rcutij; - int *element; // index on [0,nelements) - int nmax; + double **rij; // short rij list + int *inside; // short neighbor list + double *wj; // short weight list + double *rcutij; // short cutoff list + double *rinnerij; // short inner cutoff list + double *drinnerij;// short inner range list + int *element; // short element list [0,nelements) + int nmax; // allocated size of short lists void grow_rij(int); @@ -90,7 +92,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; @@ -107,11 +109,12 @@ class SNA : protected Pointers { void print_clebsch_gordan(); void init_rootpqarray(); void zero_uarraytot(int); - void add_uarraytot(double, double, double, int, int); + void add_uarraytot(double, double, double, int, int, double, double); 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, double, double); // Sets the style for the switching function // 0 = none From 9e78818637f7c504544775acb6e1141a1acc1fb9 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 11:26:06 -0700 Subject: [PATCH 03/16] Minor tweaks --- src/ML-SNAP/pair_snap.cpp | 2 +- src/ML-SNAP/sna.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index e263197eb2..9d6ebad46d 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -179,7 +179,7 @@ void PairSNAP::compute(int eflag, int vflag) else snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], snaptr->rcutij[jj], jj, 0, - // snaptr->rinnerij[jj], snaptr->drinnerij[jj],jj); +// snaptr->rinnerij[jj], snaptr->drinnerij[jj],jj); 0.0, 1.0); snaptr->compute_deidrj(fij); diff --git a/src/ML-SNAP/sna.cpp b/src/ML-SNAP/sna.cpp index 58a1597a4a..2f0d994b92 100644 --- a/src/ML-SNAP/sna.cpp +++ b/src/ML-SNAP/sna.cpp @@ -376,7 +376,8 @@ void SNA::compute_ui(int jnum, int ielem) if (chem_flag) add_uarraytot(r, wj[j], rcutij[j], j, element[j], rinnerij[j], drinnerij[j]); else - add_uarraytot(r, wj[j], rcutij[j], j, 0, rinnerij[j], drinnerij[j]); + // add_uarraytot(r, wj[j], rcutij[j], j, 0, rinnerij[j], drinnerij[j]); + add_uarraytot(r, wj[j], rcutij[j], j, 0, 0.0, 1.0); } } From 3f9cc03fc0c4b6b6a13d8071f0742fcf1218188d Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 11:39:06 -0700 Subject: [PATCH 04/16] Minor tweaks --- src/ML-SNAP/pair_snap.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index 9d6ebad46d..57abbeeb28 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -151,7 +151,7 @@ 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 (chemflag) snaptr->element[ninside] = jelem; ninside++; } } @@ -328,7 +328,7 @@ 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 (chemflag) snaptr->element[ninside] = jelem; ninside++; } } From 8323c4fe4520de995396198cbd5a4b83e3c5badc Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 14:14:45 -0700 Subject: [PATCH 05/16] More progress --- src/ML-SNAP/pair_snap.cpp | 53 +++++++++++++++++++++++++++++++++++---- src/ML-SNAP/pair_snap.h | 4 ++- src/ML-SNAP/sna.cpp | 5 ++-- src/ML-SNAP/sna.h | 29 ++++++++++++++------- 4 files changed, 73 insertions(+), 18 deletions(-) diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index 57abbeeb28..9896423ccf 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; @@ -62,6 +64,8 @@ PairSNAP::~PairSNAP() memory->destroy(radelem); memory->destroy(wjelem); memory->destroy(coeffelem); + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); memory->destroy(beta); memory->destroy(bispectrum); @@ -151,6 +155,12 @@ void PairSNAP::compute(int eflag, int vflag) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); + snaptr->rinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + // snaptr->rinnerij[ninside] = 0.0; + // snaptr->drinnerij[ninside] = 1.5; + } if (chemflag) snaptr->element[ninside] = jelem; ninside++; } @@ -179,8 +189,7 @@ void PairSNAP::compute(int eflag, int vflag) else snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], snaptr->rcutij[jj], jj, 0, -// snaptr->rinnerij[jj], snaptr->drinnerij[jj],jj); - 0.0, 1.0); + snaptr->rinnerij[jj], snaptr->drinnerij[jj]); snaptr->compute_deidrj(fij); @@ -328,6 +337,12 @@ void PairSNAP::compute_bispectrum() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; + if (switchinnerflag) { + snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); + snaptr->rinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + // snaptr->rinnerij[ninside] = 0.0; + // snaptr->drinnerij[ninside] = 1.5; + } if (chemflag) snaptr->element[ninside] = jelem; ninside++; } @@ -512,9 +527,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 @@ -630,10 +649,15 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) chemflag = 0; bnormflag = 0; wselfallflag = 0; - switchinnerflag = 1; + 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; @@ -667,7 +691,7 @@ 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]; @@ -698,11 +722,25 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) 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 + else if (keywd == "rinner") { + for (int ielem = 0; ielem < nelements; ielem++) { + rinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + keyval = strtok(nullptr,"' \t\n\r\f"); + } + rinnerflag = 1; + } else if (keywd == "drinner") { + for (int ielem = 0; ielem < nelements; ielem++) { + drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + keyval = strtok(nullptr,"' \t\n\r\f"); + } + drinnerflag = 1; + } else error->all(FLERR,"Unknown parameter '{}' in SNAP parameter file", keywd); } @@ -712,6 +750,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 5cab36d1c4..77ceeabd4f 100644 --- a/src/ML-SNAP/pair_snap.h +++ b/src/ML-SNAP/pair_snap.h @@ -59,7 +59,9 @@ class PairSNAP : public Pair { double **scale; // for thermodynamic integration int twojmax, switchflag, bzeroflag, bnormflag; int chemflag, wselfallflag; - int switchinnerflag; + 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 2f0d994b92..39d50bcab7 100644 --- a/src/ML-SNAP/sna.cpp +++ b/src/ML-SNAP/sna.cpp @@ -178,7 +178,7 @@ SNA::~SNA() memory->destroy(rinnerij); memory->destroy(drinnerij); } -if (chem_flag) memory->destroy(element); + if (chem_flag) memory->destroy(element); memory->destroy(ulist_r_ij); memory->destroy(ulist_i_ij); delete[] idxz; @@ -376,8 +376,7 @@ void SNA::compute_ui(int jnum, int ielem) if (chem_flag) add_uarraytot(r, wj[j], rcutij[j], j, element[j], rinnerij[j], drinnerij[j]); else - // add_uarraytot(r, wj[j], rcutij[j], j, 0, rinnerij[j], drinnerij[j]); - add_uarraytot(r, wj[j], rcutij[j], j, 0, 0.0, 1.0); + add_uarraytot(r, wj[j], rcutij[j], j, 0, rinnerij[j], drinnerij[j]); } } diff --git a/src/ML-SNAP/sna.h b/src/ML-SNAP/sna.h index 8c631abccd..bad0db1a7c 100644 --- a/src/ML-SNAP/sna.h +++ b/src/ML-SNAP/sna.h @@ -61,23 +61,31 @@ class SNA : protected Pointers { double compute_sfac_inner(double, double, double); double compute_dsfac_inner(double, double, double); - + + // public bispectrum data + + int twojmax; 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 + 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) - int nmax; // allocated size of short lists - - void grow_rij(int); - - int twojmax; - double *ylist_r, *ylist_i; - int idxcg_max, idxu_max, idxz_max, idxb_max; private: double rmin0, rfac0; @@ -103,6 +111,9 @@ 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(); From 877b764a362c0369bd575c01ddb0d38d21c5f155 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 14:38:51 -0700 Subject: [PATCH 06/16] More --- src/ML-SNAP/sna.cpp | 20 ++++++++++---------- src/ML-SNAP/sna.h | 2 +- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/ML-SNAP/sna.cpp b/src/ML-SNAP/sna.cpp index 39d50bcab7..8527f42b3e 100644 --- a/src/ML-SNAP/sna.cpp +++ b/src/ML-SNAP/sna.cpp @@ -373,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], rinnerij[j], drinnerij[j]); - else - add_uarraytot(r, wj[j], rcutij[j], j, 0, rinnerij[j], drinnerij[j]); + add_uarraytot(r, j); } } @@ -1014,16 +1011,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, - double rinner, double drinner) +void SNA::add_uarraytot(double r, int jj) { double sfac; - - sfac = compute_sfac(r, rcut); - sfac *= wj; + int jelem; + + sfac = compute_sfac(r, rcutij[jj]); + sfac *= wj[jj]; if (switch_inner_flag) - sfac *= compute_sfac_inner(r, rinner, drinner); + 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]; diff --git a/src/ML-SNAP/sna.h b/src/ML-SNAP/sna.h index bad0db1a7c..e2fa06d803 100644 --- a/src/ML-SNAP/sna.h +++ b/src/ML-SNAP/sna.h @@ -120,7 +120,7 @@ class SNA : protected Pointers { void print_clebsch_gordan(); void init_rootpqarray(); void zero_uarraytot(int); - void add_uarraytot(double, double, double, int, int, double, double); + void add_uarraytot(double, int); void compute_uarray(double, double, double, double, double, int); double deltacg(int, int, int); void compute_ncoeff(); From 4d132821208acbbf23c06a6f35c02644f34359a7 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 15:07:18 -0700 Subject: [PATCH 07/16] More --- src/ML-IAP/mliap_descriptor_snap.cpp | 29 +++------------------------- src/ML-SNAP/compute_snad_atom.cpp | 4 +--- src/ML-SNAP/compute_snap.cpp | 4 +--- src/ML-SNAP/compute_snav_atom.cpp | 4 +--- src/ML-SNAP/pair_snap.cpp | 9 +-------- src/ML-SNAP/sna.cpp | 23 +++++++++++----------- src/ML-SNAP/sna.h | 4 ++-- 7 files changed, 21 insertions(+), 56 deletions(-) diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index 6571637681..9684a056e9 100644 --- a/src/ML-IAP/mliap_descriptor_snap.cpp +++ b/src/ML-IAP/mliap_descriptor_snap.cpp @@ -161,14 +161,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], - snaptr->rinnerij[jj], snaptr->drinnerij[jj]); - else - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj], jj, 0, - snaptr->rinnerij[jj], snaptr->drinnerij[jj]); + snaptr->compute_duidrj(jj); snaptr->compute_deidrj(fij); @@ -237,15 +230,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], - snaptr->rinnerij[jj], snaptr->drinnerij[jj]); - else - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj], jj, 0, - snaptr->rinnerij[jj], snaptr->drinnerij[jj]); - + snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate gamma_lk*dB_k/dRi, -gamma_lk**dB_k/dRj @@ -312,15 +297,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], - snaptr->rinnerij[jj], snaptr->drinnerij[jj]); - else - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj], jj, 0, - snaptr->rinnerij[jj], snaptr->drinnerij[jj]); - + snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate dB_k^i/dRi, dB_k^i/dRj diff --git a/src/ML-SNAP/compute_snad_atom.cpp b/src/ML-SNAP/compute_snad_atom.cpp index a94555faa2..acee813dae 100644 --- a/src/ML-SNAP/compute_snad_atom.cpp +++ b/src/ML-SNAP/compute_snad_atom.cpp @@ -301,9 +301,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->rinnerij[jj], snaptr->drinnerij[jj]); + snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate -dBi/dRi, -dBi/dRj diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp index 6cf7a38235..f371f3ff55 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -355,9 +355,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->rinnerij[jj], snaptr->drinnerij[jj]); + snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate dBi/dRi, -dBi/dRj diff --git a/src/ML-SNAP/compute_snav_atom.cpp b/src/ML-SNAP/compute_snav_atom.cpp index dd49a90b90..6b53b344e8 100644 --- a/src/ML-SNAP/compute_snav_atom.cpp +++ b/src/ML-SNAP/compute_snav_atom.cpp @@ -295,9 +295,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->rinnerij[jj], snaptr->drinnerij[jj]); + snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate -dBi/dRi*Ri, -dBi/dRj*Rj diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index 9896423ccf..1d4db2faae 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -182,14 +182,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], - snaptr->rinnerij[jj], snaptr->drinnerij[jj]); - else - snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], - snaptr->rcutij[jj], jj, 0, - snaptr->rinnerij[jj], snaptr->drinnerij[jj]); + snaptr->compute_duidrj(jj); snaptr->compute_deidrj(fij); diff --git a/src/ML-SNAP/sna.cpp b/src/ML-SNAP/sna.cpp index 8527f42b3e..9488f83112 100644 --- a/src/ML-SNAP/sna.cpp +++ b/src/ML-SNAP/sna.cpp @@ -963,15 +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, double rinner, double drinner) +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); @@ -981,8 +981,10 @@ void SNA::compute_duidrj(double* rij, double wj, double rcut, int jj, 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, rinner, drinner); + if (chem_flag) elem_duarray = element[jj]; + else elem_duarray = 0; + + compute_duarray(x, y, z, z0, r, dz0dr, wj[jj], rcut, jj); } /* ---------------------------------------------------------------------- */ @@ -1137,8 +1139,7 @@ void SNA::compute_uarray(double x, double y, double z, void SNA::compute_duarray(double x, double y, double z, double z0, double r, double dz0dr, - double wj, double rcut, int jj, - double rinner, double drinner) + double wj, double rcut, int jj) { double r0inv; double a_r, a_i, b_r, b_i; @@ -1269,8 +1270,8 @@ void SNA::compute_duarray(double x, double y, double z, double dsfac = compute_dsfac(r, rcut); if (switch_inner_flag) { - double sfac_inner = compute_sfac_inner(r, rinner, drinner); - dsfac = dsfac*sfac_inner + sfac*compute_dsfac_inner(r, rinner, drinner); + 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; } diff --git a/src/ML-SNAP/sna.h b/src/ML-SNAP/sna.h index e2fa06d803..126410fa00 100644 --- a/src/ML-SNAP/sna.h +++ b/src/ML-SNAP/sna.h @@ -53,7 +53,7 @@ class SNA : protected Pointers { // functions for derivatives - void compute_duidrj(double *, double, double, int, int, double, double); + void compute_duidrj(int); void compute_dbidrj(); void compute_deidrj(double *); double compute_sfac(double, double); @@ -125,7 +125,7 @@ class SNA : protected Pointers { double deltacg(int, int, int); void compute_ncoeff(); void compute_duarray(double, double, double, double, double, double, - double, double, int, double, double); + double, double, int); // Sets the style for the switching function // 0 = none From 5216a543ae2b14abefb4f2461612b919c5100444 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 15:43:48 -0700 Subject: [PATCH 08/16] More --- src/ML-IAP/mliap_descriptor_snap.cpp | 19 +++++++++++-------- src/ML-SNAP/compute_sna_atom.cpp | 2 +- src/ML-SNAP/compute_snad_atom.cpp | 2 +- src/ML-SNAP/compute_snap.cpp | 2 +- src/ML-SNAP/compute_snav_atom.cpp | 2 +- 5 files changed, 15 insertions(+), 12 deletions(-) diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index 9684a056e9..aa2b2d34af 100644 --- a/src/ML-IAP/mliap_descriptor_snap.cpp +++ b/src/ML-IAP/mliap_descriptor_snap.cpp @@ -88,7 +88,7 @@ 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 (chemflag) snaptr->element[ninside] = jelem; ninside++; ij++; } @@ -141,7 +141,7 @@ 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 (chemflag) snaptr->element[ninside] = jelem; ninside++; ij++; } @@ -211,7 +211,7 @@ 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 (chemflag) snaptr->element[ninside] = jelem; ninside++; ij++; } @@ -279,7 +279,7 @@ 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 (chemflag) snaptr->element[ninside] = jelem; ninside++; ij++; } @@ -392,9 +392,6 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) 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); - // check for keywords with one value per element if (strcmp(keywd,"elems") == 0 || @@ -404,6 +401,9 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) if (nelementsflag == 0 || nwords != nelements+1) error->all(FLERR,"Incorrect SNAP parameter file"); + if (comm->me == 0) + utils::logmesg(lmp,"SNAP keyword {} {} ... \n", keywd, keyval); + if (strcmp(keywd,"elems") == 0) { for (int ielem = 0; ielem < nelements; ielem++) { elements[ielem] = utils::strdup(keyval); @@ -426,11 +426,14 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) } else { - // all other keywords take one value + // 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 (strcmp(keywd,"nelems") == 0) { nelements = atoi(keyval); elements = new char*[nelements]; diff --git a/src/ML-SNAP/compute_sna_atom.cpp b/src/ML-SNAP/compute_sna_atom.cpp index 5dd337286b..421a317364 100644 --- a/src/ML-SNAP/compute_sna_atom.cpp +++ b/src/ML-SNAP/compute_sna_atom.cpp @@ -270,7 +270,7 @@ 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 (chemflag) snaptr->element[ninside] = jelem; ninside++; } } diff --git a/src/ML-SNAP/compute_snad_atom.cpp b/src/ML-SNAP/compute_snad_atom.cpp index acee813dae..2f2ad6f970 100644 --- a/src/ML-SNAP/compute_snad_atom.cpp +++ b/src/ML-SNAP/compute_snad_atom.cpp @@ -288,7 +288,7 @@ 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 (chemflag) snaptr->element[ninside] = jelem; ninside++; } } diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp index f371f3ff55..6d4649ed4f 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -344,7 +344,7 @@ 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 (chemflag) snaptr->element[ninside] = jelem; ninside++; } } diff --git a/src/ML-SNAP/compute_snav_atom.cpp b/src/ML-SNAP/compute_snav_atom.cpp index 6b53b344e8..487838a770 100644 --- a/src/ML-SNAP/compute_snav_atom.cpp +++ b/src/ML-SNAP/compute_snav_atom.cpp @@ -281,7 +281,7 @@ 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 (chemflag) snaptr->element[ninside] = jelem; ninside++; } } From f691805062c9b13ae7897503b51bbc060a90c4b5 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 16:41:22 -0700 Subject: [PATCH 09/16] More --- src/ML-IAP/mliap_descriptor_snap.cpp | 2 +- src/ML-SNAP/compute_sna_atom.cpp | 2 +- src/ML-SNAP/compute_snad_atom.cpp | 2 +- src/ML-SNAP/compute_snap.cpp | 2 +- src/ML-SNAP/compute_snav_atom.cpp | 2 +- src/ML-SNAP/pair_snap.cpp | 112 ++++++++++++++++----------- 6 files changed, 73 insertions(+), 49 deletions(-) diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index aa2b2d34af..5dbb461c09 100644 --- a/src/ML-IAP/mliap_descriptor_snap.cpp +++ b/src/ML-IAP/mliap_descriptor_snap.cpp @@ -345,7 +345,7 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) chemflag = 0; bnormflag = 0; wselfallflag = 0; - switchinnerflag = 1; + switchinnerflag = 0; for (int i = 0; i < nelements; i++) delete[] elements[i]; delete[] elements; diff --git a/src/ML-SNAP/compute_sna_atom.cpp b/src/ML-SNAP/compute_sna_atom.cpp index 421a317364..4a3094a9fd 100644 --- a/src/ML-SNAP/compute_sna_atom.cpp +++ b/src/ML-SNAP/compute_sna_atom.cpp @@ -54,7 +54,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; - switchinnerflag = 1; + switchinnerflag = 0; nelements = 1; // offset by 1 to match up with types diff --git a/src/ML-SNAP/compute_snad_atom.cpp b/src/ML-SNAP/compute_snad_atom.cpp index 2f2ad6f970..1d92065f2a 100644 --- a/src/ML-SNAP/compute_snad_atom.cpp +++ b/src/ML-SNAP/compute_snad_atom.cpp @@ -54,7 +54,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; - switchinnerflag = 1; + switchinnerflag = 0; nelements = 1; // process required arguments diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp index 6d4649ed4f..da870a1728 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -61,7 +61,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; - switchinnerflag = 1; + switchinnerflag = 0; nelements = 1; // process required arguments diff --git a/src/ML-SNAP/compute_snav_atom.cpp b/src/ML-SNAP/compute_snav_atom.cpp index 487838a770..b1d2d2b51b 100644 --- a/src/ML-SNAP/compute_snav_atom.cpp +++ b/src/ML-SNAP/compute_snav_atom.cpp @@ -53,7 +53,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : chemflag = 0; bnormflag = 0; wselfallflag = 0; - switchinnerflag = 1; + switchinnerflag = 0; nelements = 1; // process required arguments diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index 1d4db2faae..8544b728e0 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -674,7 +674,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 { @@ -684,57 +685,80 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) } if (words.size() == 0) continue; + 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 == "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 if (keywd == "rinner") { - for (int ielem = 0; ielem < nelements; ielem++) { - rinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); - keyval = strtok(nullptr,"' \t\n\r\f"); + 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); + + if (keywd == "rinner") { + for (int ielem = 0; ielem < nelements; ielem++) { + rinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + keyval = strtok(nullptr,"' \t\n\r\f"); + } + rinnerflag = 1; + } else if (keywd == "drinner") { + printf("drinner nelements = %d %s\n",nelements,keyval.c_str()); + for (int ielem = 0; ielem < nelements; ielem++) { + drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + keyval = strtok(nullptr,"' \t\n\r\f"); + printf("drinner ielem = %d drinnerelem[ielem] = %g\n",ielem,drinnerelem[ielem]); + } + drinnerflag = 1; } - rinnerflag = 1; - } else if (keywd == "drinner") { - for (int ielem = 0; ielem < nelements; ielem++) { - drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); - keyval = strtok(nullptr,"' \t\n\r\f"); - } - drinnerflag = 1; - } else - error->all(FLERR,"Unknown parameter '{}' in SNAP parameter file", keywd); + + } 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) From b9fd86e7ca041b91484bded791758298e6c78303 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 17:23:38 -0700 Subject: [PATCH 10/16] switchinnerflag working for single-element two types SNAP potential --- src/ML-SNAP/pair_snap.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index 8544b728e0..0430b36617 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -157,7 +157,7 @@ void PairSNAP::compute(int eflag, int vflag) snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; if (switchinnerflag) { snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); - snaptr->rinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); // snaptr->rinnerij[ninside] = 0.0; // snaptr->drinnerij[ninside] = 1.5; } @@ -332,7 +332,7 @@ void PairSNAP::compute_bispectrum() snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; if (switchinnerflag) { snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); - snaptr->rinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); // snaptr->rinnerij[ninside] = 0.0; // snaptr->drinnerij[ninside] = 1.5; } @@ -702,18 +702,20 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) 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++) { rinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); - keyval = strtok(nullptr,"' \t\n\r\f"); + iword++; } rinnerflag = 1; } else if (keywd == "drinner") { - printf("drinner nelements = %d %s\n",nelements,keyval.c_str()); + keyval = words[iword]; for (int ielem = 0; ielem < nelements; ielem++) { drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); - keyval = strtok(nullptr,"' \t\n\r\f"); - printf("drinner ielem = %d drinnerelem[ielem] = %g\n",ielem,drinnerelem[ielem]); + iword++; } drinnerflag = 1; } From 51034c00c089daff277c81f83db8cd14f33a2a15 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 20:40:18 -0700 Subject: [PATCH 11/16] Completed documentation --- doc/src/compute_sna_atom.rst | 26 +++++++- doc/src/pair_snap.rst | 14 +++- doc/utils/sphinx-config/false_positives.txt | 2 + src/ML-IAP/mliap_descriptor_snap.cpp | 74 +++++++++++++++++---- src/ML-IAP/mliap_descriptor_snap.h | 3 + src/ML-SNAP/compute_sna_atom.cpp | 31 +++++++-- src/ML-SNAP/compute_sna_atom.h | 2 + src/ML-SNAP/compute_snad_atom.cpp | 28 +++++++- src/ML-SNAP/compute_snad_atom.h | 2 + src/ML-SNAP/compute_snap.cpp | 26 +++++++- src/ML-SNAP/compute_snap.h | 2 + src/ML-SNAP/compute_snav_atom.cpp | 29 ++++++-- src/ML-SNAP/compute_snav_atom.h | 2 + src/ML-SNAP/pair_snap.cpp | 29 ++++---- 14 files changed, 229 insertions(+), 41 deletions(-) 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 e033b9088f..d996d46346 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1927,6 +1927,7 @@ mbt MBytes mc McLachlan +mcmoves md mdf mdi @@ -2227,6 +2228,7 @@ Nelement Nelements nelems nemd +netapp netcdf netstat Nettleton diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index 5dbb461c09..6de7c0c3ca 100644 --- a/src/ML-IAP/mliap_descriptor_snap.cpp +++ b/src/ML-IAP/mliap_descriptor_snap.cpp @@ -42,6 +42,9 @@ MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename): radelem = nullptr; wjelem = nullptr; snaptr = nullptr; + rinnerelem = nullptr; + drinnerelem = nullptr; + read_paramfile(paramfilename); snaptr = new SNA(lmp, rfac0, twojmax, @@ -59,6 +62,11 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP() memory->destroy(radelem); memory->destroy(wjelem); delete snaptr; + + if (switchinnerflag) { + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); + } } /* ---------------------------------------------------------------------- @@ -88,6 +96,10 @@ void MLIAPDescriptorSNAP::compute_descriptors(class MLIAPData* data) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][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++; ij++; @@ -141,6 +153,10 @@ void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][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++; ij++; @@ -211,6 +227,10 @@ void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][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++; ij++; @@ -279,6 +299,10 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][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++; ij++; @@ -347,6 +371,11 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) 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; memory->destroy(radelem); @@ -396,7 +425,9 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) if (strcmp(keywd,"elems") == 0 || strcmp(keywd,"radelems") == 0 || - strcmp(keywd,"welems") == 0) { + strcmp(keywd,"welems") == 0 || + strcmp(keywd,"rinnerelems") == 0 || + strcmp(keywd,"drinnerelems") == 0) { if (nelementsflag == 0 || nwords != nelements+1) error->all(FLERR,"Incorrect SNAP parameter file"); @@ -422,6 +453,18 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) keyval = strtok(nullptr,"' \t\n\r\f"); } wjelemflag = 1; + } else if (strcmp(keywd,"rinnerelems") == 0) { + for (int ielem = 0; ielem < nelements; ielem++) { + rinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + keyval = strtok(nullptr,"' \t\n\r\f"); + } + rinnerflag = 1; + } else if (strcmp(keywd,"drinnerelems") == 0) { + for (int ielem = 0; ielem < nelements; ielem++) { + drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + keyval = strtok(nullptr,"' \t\n\r\f"); + } + rinnerflag = 1; } } else { @@ -435,34 +478,35 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) utils::logmesg(lmp,"SNAP keyword {} {} \n", keywd, keyval); if (strcmp(keywd,"nelems") == 0) { - nelements = atoi(keyval); + 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"); nelementsflag = 1; } else if (strcmp(keywd,"rcutfac") == 0) { - rcutfac = atof(keyval); + rcutfac = utils::numeric(FLERR,keyval,false,lmp); rcutfacflag = 1; } else if (strcmp(keywd,"twojmax") == 0) { - twojmax = atoi(keyval); + twojmax = utils::inumeric(FLERR,keyval,false,lmp); twojmaxflag = 1; } else if (strcmp(keywd,"rfac0") == 0) - rfac0 = atof(keyval); + rfac0 = utils::numeric(FLERR,keyval,false,lmp); else if (strcmp(keywd,"rmin0") == 0) - rmin0 = atof(keyval); + rmin0 = utils::numeric(FLERR,keyval,false,lmp); else if (strcmp(keywd,"switchflag") == 0) - switchflag = atoi(keyval); + switchflag = utils::inumeric(FLERR,keyval,false,lmp); else if (strcmp(keywd,"bzeroflag") == 0) - bzeroflag = atoi(keyval); + bzeroflag = utils::inumeric(FLERR,keyval,false,lmp); else if (strcmp(keywd,"chemflag") == 0) - chemflag = atoi(keyval); + chemflag = utils::inumeric(FLERR,keyval,false,lmp); else if (strcmp(keywd,"bnormflag") == 0) - bnormflag = atoi(keyval); + bnormflag = utils::inumeric(FLERR,keyval,false,lmp); else if (strcmp(keywd,"wselfallflag") == 0) - wselfallflag = atoi(keyval); + 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"); - } } @@ -470,6 +514,12 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) !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; diff --git a/src/ML-IAP/mliap_descriptor_snap.h b/src/ML-IAP/mliap_descriptor_snap.h index 8c8a8f16f4..2136eefc5e 100644 --- a/src/ML-IAP/mliap_descriptor_snap.h +++ b/src/ML-IAP/mliap_descriptor_snap.h @@ -41,6 +41,9 @@ class MLIAPDescriptorSNAP : public MLIAPDescriptor { int chemflag, bnormflag, wselfallflag; int switchinnerflag; double rfac0, rmin0; + + double* rinnerelem; + double* drinnerelem; }; } // namespace LAMMPS_NS diff --git a/src/ML-SNAP/compute_sna_atom.cpp b/src/ML-SNAP/compute_sna_atom.cpp index 4a3094a9fd..8a99c0e14d 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; @@ -134,7 +133,20 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute sna/atom command"); wselfallflag = atoi(arg[iarg+1]); iarg += 2; - } else error->all(FLERR,"Illegal compute sna/atom command"); + } 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, @@ -160,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); + } } /* ---------------------------------------------------------------------- */ @@ -270,6 +289,10 @@ void ComputeSNAAtom::compute_peratom() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; + 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 e741e02bc2..5e62fe47d9 100644 --- a/src/ML-SNAP/compute_sna_atom.h +++ b/src/ML-SNAP/compute_sna_atom.h @@ -45,6 +45,8 @@ class ComputeSNAAtom : public Compute { 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 1d92065f2a..0fb24c516a 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; @@ -132,6 +130,19 @@ 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"); } @@ -162,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); + } } /* ---------------------------------------------------------------------- */ @@ -288,6 +306,10 @@ void ComputeSNADAtom::compute_peratom() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; + 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_snad_atom.h b/src/ML-SNAP/compute_snad_atom.h index 7727a4653b..c03458446f 100644 --- a/src/ML-SNAP/compute_snad_atom.h +++ b/src/ML-SNAP/compute_snad_atom.h @@ -47,6 +47,8 @@ class ComputeSNADAtom : public Compute { 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 da870a1728..d9bc4afb96 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; @@ -144,6 +142,19 @@ 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"); } @@ -185,6 +196,11 @@ ComputeSnap::~ComputeSnap() delete snaptr; if (chemflag) memory->destroy(map); + + if (switchinnerflag) { + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); + } } /* ---------------------------------------------------------------------- */ @@ -344,6 +360,10 @@ void ComputeSnap::compute_array() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; + 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_snap.h b/src/ML-SNAP/compute_snap.h index 9d928a8166..fa31de747d 100644 --- a/src/ML-SNAP/compute_snap.h +++ b/src/ML-SNAP/compute_snap.h @@ -47,6 +47,8 @@ class ComputeSnap : public Compute { 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 b1d2d2b51b..a6cb751eda 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; @@ -127,6 +125,19 @@ 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"); } @@ -155,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); + } } /* ---------------------------------------------------------------------- */ @@ -281,6 +298,10 @@ void ComputeSNAVAtom::compute_peratom() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; + 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_snav_atom.h b/src/ML-SNAP/compute_snav_atom.h index ea5ea50086..be442f3a49 100644 --- a/src/ML-SNAP/compute_snav_atom.h +++ b/src/ML-SNAP/compute_snav_atom.h @@ -47,6 +47,8 @@ class ComputeSNAVAtom : public Compute { 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 0430b36617..ba8a498411 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -64,9 +64,12 @@ PairSNAP::~PairSNAP() memory->destroy(radelem); memory->destroy(wjelem); memory->destroy(coeffelem); - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); + if (switchinnerflag) { + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); + } + memory->destroy(beta); memory->destroy(bispectrum); @@ -158,8 +161,6 @@ void PairSNAP::compute(int eflag, int vflag) if (switchinnerflag) { snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); - // snaptr->rinnerij[ninside] = 0.0; - // snaptr->drinnerij[ninside] = 1.5; } if (chemflag) snaptr->element[ninside] = jelem; ninside++; @@ -333,8 +334,6 @@ void PairSNAP::compute_bispectrum() if (switchinnerflag) { snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); - // snaptr->rinnerij[ninside] = 0.0; - // snaptr->drinnerij[ninside] = 1.5; } if (chemflag) snaptr->element[ninside] = jelem; ninside++; @@ -520,14 +519,18 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) memory->destroy(radelem); memory->destroy(wjelem); memory->destroy(coeffelem); - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); + if (switchinnerflag) { + 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"); - + if (switchinnerflag) { + memory->create(rinnerelem,nelements,"pair:rinnerelem"); + memory->create(drinnerelem,nelements,"pair:drinnerelem"); + } + // initialize checklist for all required nelements int *elementflags = new int[nelements]; @@ -707,14 +710,14 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) if (keywd == "rinner") { keyval = words[iword]; for (int ielem = 0; ielem < nelements; ielem++) { - rinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + 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); + drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); iword++; } drinnerflag = 1; From 014345a5126c352bce86bf30fca628bc21e4dde7 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 20:55:59 -0700 Subject: [PATCH 12/16] Updated a citation --- potentials/InP_JCPA2020.mliap | 2 +- potentials/InP_JCPA2020.mliap.descriptor | 2 +- potentials/InP_JCPA2020.mliap.model | 2 +- potentials/InP_JCPA2020.snap | 2 +- potentials/InP_JCPA2020.snapcoeff | 2 +- potentials/InP_JCPA2020.snapparam | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) 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 From b1119d81a96352463767656ee19f5ef9a73311e3 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Mar 2022 21:25:44 -0700 Subject: [PATCH 13/16] Minor tweaks --- src/ML-IAP/mliap_descriptor_snap.cpp | 10 ++++++---- src/ML-SNAP/pair_snap.cpp | 13 +++++-------- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index 6de7c0c3ca..37387cc7d7 100644 --- a/src/ML-IAP/mliap_descriptor_snap.cpp +++ b/src/ML-IAP/mliap_descriptor_snap.cpp @@ -64,9 +64,9 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP() delete snaptr; if (switchinnerflag) { - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); - } + memory->destroy(rinnerelem); + memory->destroy(drinnerelem); + } } /* ---------------------------------------------------------------------- @@ -482,6 +482,8 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) 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 = utils::numeric(FLERR,keyval,false,lmp); @@ -503,7 +505,7 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) bnormflag = utils::inumeric(FLERR,keyval,false,lmp); else if (strcmp(keywd,"wselfallflag") == 0) wselfallflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (keywd == "switchinnerflag") + else if (strcmp(keywd,"switchinnerflag") == 0) switchinnerflag = utils::inumeric(FLERR,keyval,false,lmp); else error->all(FLERR,"Incorrect SNAP parameter file"); diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index ba8a498411..5166c5e0e9 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -519,17 +519,13 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) memory->destroy(radelem); memory->destroy(wjelem); memory->destroy(coeffelem); - if (switchinnerflag) { - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); - } + 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"); - if (switchinnerflag) { - memory->create(rinnerelem,nelements,"pair:rinnerelem"); - memory->create(drinnerelem,nelements,"pair:drinnerelem"); - } + memory->create(rinnerelem,nelements,"pair:rinnerelem"); + memory->create(drinnerelem,nelements,"pair:drinnerelem"); // initialize checklist for all required nelements @@ -710,6 +706,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) 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++; } From 1d1dea56d8cbb533700d629398d1cbad93d20a0e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 14 Mar 2022 16:52:17 -0400 Subject: [PATCH 14/16] whitespace and programming style fixes --- src/ML-IAP/mliap_descriptor_snap.cpp | 103 ++++++++++++--------------- src/ML-SNAP/compute_sna_atom.cpp | 12 ++-- src/ML-SNAP/compute_snad_atom.cpp | 10 +-- src/ML-SNAP/compute_snap.cpp | 14 ++-- src/ML-SNAP/compute_snav_atom.cpp | 10 +-- src/ML-SNAP/pair_snap.cpp | 88 +++++++++++------------ src/ML-SNAP/sna.cpp | 10 +-- src/ML-SNAP/sna.h | 8 +-- 8 files changed, 123 insertions(+), 132 deletions(-) diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index 37387cc7d7..a914e523e2 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 @@ -36,8 +35,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename): - MLIAPDescriptor(lmp) +MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *_lmp, char *paramfilename) : MLIAPDescriptor(_lmp) { radelem = nullptr; wjelem = nullptr; @@ -47,10 +45,8 @@ MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename): read_paramfile(paramfilename); - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); ndescriptors = snaptr->ncoeff; } @@ -73,7 +69,7 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP() 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++) { @@ -97,8 +93,8 @@ void MLIAPDescriptorSNAP::compute_descriptors(class MLIAPData* data) snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + 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++; @@ -119,14 +115,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; @@ -154,8 +149,8 @@ void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data) snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + 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++; @@ -191,18 +186,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++) { @@ -228,8 +221,8 @@ void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data) snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + 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++; @@ -258,24 +251,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++) { @@ -300,8 +291,8 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data) snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + 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++; @@ -334,7 +325,6 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data) ij++; } } - } /* ---------------------------------------------------------------------- @@ -386,28 +376,29 @@ 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 @@ -512,28 +503,28 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) } } - 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"); - 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; } } } @@ -545,7 +536,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-SNAP/compute_sna_atom.cpp b/src/ML-SNAP/compute_sna_atom.cpp index 8a99c0e14d..0dc2de4656 100644 --- a/src/ML-SNAP/compute_sna_atom.cpp +++ b/src/ML-SNAP/compute_sna_atom.cpp @@ -146,13 +146,13 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : 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"); + } else error->all(FLERR,"Illegal compute sna/atom command"); } snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; size_peratom_cols = ncoeff; @@ -289,10 +289,10 @@ void ComputeSNAAtom::compute_peratom() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; - if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); - } + 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_snad_atom.cpp b/src/ML-SNAP/compute_snad_atom.cpp index 0fb24c516a..25f4b430aa 100644 --- a/src/ML-SNAP/compute_snad_atom.cpp +++ b/src/ML-SNAP/compute_snad_atom.cpp @@ -149,7 +149,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -306,10 +306,10 @@ void ComputeSNADAtom::compute_peratom() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; - if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); - } + 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_snap.cpp b/src/ML-SNAP/compute_snap.cpp index d9bc4afb96..af8c52c2ee 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -150,10 +150,10 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : 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); + 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); + drinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); iarg += ntypes; } else error->all(FLERR,"Illegal compute snap command"); } @@ -161,7 +161,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -360,10 +360,10 @@ void ComputeSnap::compute_array() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; - if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); - } + 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_snav_atom.cpp b/src/ML-SNAP/compute_snav_atom.cpp index a6cb751eda..9834dc2a0c 100644 --- a/src/ML-SNAP/compute_snav_atom.cpp +++ b/src/ML-SNAP/compute_snav_atom.cpp @@ -144,7 +144,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + nelements, switchinnerflag); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -298,10 +298,10 @@ void ComputeSNAVAtom::compute_peratom() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; - if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); - } + 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/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index 5166c5e0e9..e9c6958ceb 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -69,7 +69,7 @@ PairSNAP::~PairSNAP() memory->destroy(rinnerelem); memory->destroy(drinnerelem); } - + memory->destroy(beta); memory->destroy(bispectrum); @@ -158,11 +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; - 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; + 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++; } } @@ -331,11 +331,11 @@ void PairSNAP::compute_bispectrum() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; - 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; + 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++; } } @@ -414,7 +414,7 @@ void PairSNAP::coeff(int narg, char **arg) snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + nelements, switchinnerflag); if (ncoeff != snaptr->ncoeff) { if (comm->me == 0) @@ -526,7 +526,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) 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 int *elementflags = new int[nelements]; @@ -649,7 +649,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) int rinnerflag = 0; int drinnerflag = 0; - + // open SNAP parameter file on proc 0 FILE *fpparam; @@ -699,23 +699,23 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) error->all(FLERR,"Incorrect SNAP parameter file"); if (comm->me == 0) - utils::logmesg(lmp,"SNAP keyword {} {} ... \n", keywd, keyval); + utils::logmesg(lmp,"SNAP keyword {} {} ... \n", keywd, keyval); int iword = 1; if (keywd == "rinner") { - keyval = words[iword]; + 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++; + 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]; + keyval = words[iword]; for (int ielem = 0; ielem < nelements; ielem++) { - drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); - iword++; + drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + iword++; } drinnerflag = 1; } @@ -723,43 +723,43 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) } 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); - + utils::logmesg(lmp,"SNAP keyword {} {}\n",keywd,keyval); + if (keywd == "rcutfac") { - rcutfac = utils::numeric(FLERR,keyval,false,lmp); - rcutfacflag = 1; + rcutfac = utils::numeric(FLERR,keyval,false,lmp); + rcutfacflag = 1; } else if (keywd == "twojmax") { - twojmax = utils::inumeric(FLERR,keyval,false,lmp); - twojmaxflag = 1; + twojmax = utils::inumeric(FLERR,keyval,false,lmp); + twojmaxflag = 1; } else if (keywd == "rfac0") - rfac0 = utils::numeric(FLERR,keyval,false,lmp); + rfac0 = utils::numeric(FLERR,keyval,false,lmp); else if (keywd == "rmin0") - rmin0 = utils::numeric(FLERR,keyval,false,lmp); + rmin0 = utils::numeric(FLERR,keyval,false,lmp); else if (keywd == "switchflag") - switchflag = utils::inumeric(FLERR,keyval,false,lmp); + switchflag = utils::inumeric(FLERR,keyval,false,lmp); else if (keywd == "bzeroflag") - bzeroflag = utils::inumeric(FLERR,keyval,false,lmp); + bzeroflag = utils::inumeric(FLERR,keyval,false,lmp); else if (keywd == "quadraticflag") - quadraticflag = utils::inumeric(FLERR,keyval,false,lmp); + quadraticflag = utils::inumeric(FLERR,keyval,false,lmp); else if (keywd == "chemflag") - chemflag = utils::inumeric(FLERR,keyval,false,lmp); + chemflag = utils::inumeric(FLERR,keyval,false,lmp); else if (keywd == "bnormflag") - bnormflag = utils::inumeric(FLERR,keyval,false,lmp); + bnormflag = utils::inumeric(FLERR,keyval,false,lmp); else if (keywd == "wselfallflag") - wselfallflag = utils::inumeric(FLERR,keyval,false,lmp); + wselfallflag = utils::inumeric(FLERR,keyval,false,lmp); else if (keywd == "switchinnerflag") - switchinnerflag = utils::inumeric(FLERR,keyval,false,lmp); + switchinnerflag = utils::inumeric(FLERR,keyval,false,lmp); else if (keywd == "chunksize") - chunksize = utils::inumeric(FLERR,keyval,false,lmp); + chunksize = utils::inumeric(FLERR,keyval,false,lmp); else if (keywd == "parallelthresh") - parallel_thresh = utils::inumeric(FLERR,keyval,false,lmp); + parallel_thresh = utils::inumeric(FLERR,keyval,false,lmp); else - error->all(FLERR,"Unknown parameter '{}' in SNAP parameter file", keywd); + error->all(FLERR,"Unknown parameter '{}' in SNAP parameter file", keywd); } } @@ -769,10 +769,10 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) if (chemflag && nelemtmp != nelements) error->all(FLERR,"Incorrect SNAP parameter file"); - if (switchinnerflag && !(rinnerflag && drinnerflag)) + if (switchinnerflag && !(rinnerflag && drinnerflag)) error->all(FLERR,"Incorrect SNAP parameter file"); - if (!switchinnerflag && (rinnerflag || drinnerflag)) + if (!switchinnerflag && (rinnerflag || drinnerflag)) error->all(FLERR,"Incorrect SNAP parameter file"); } diff --git a/src/ML-SNAP/sna.cpp b/src/ML-SNAP/sna.cpp index 9488f83112..b33ad8a491 100644 --- a/src/ML-SNAP/sna.cpp +++ b/src/ML-SNAP/sna.cpp @@ -113,7 +113,7 @@ 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, int switch_inner_flag_in) : Pointers(lmp) + int nelements_in, int switch_inner_flag_in) : Pointers(lmp) { wself = 1.0; @@ -1017,7 +1017,7 @@ void SNA::add_uarraytot(double r, int jj) { double sfac; int jelem; - + sfac = compute_sfac(r, rcutij[jj]); sfac *= wj[jj]; @@ -1026,7 +1026,7 @@ void SNA::add_uarraytot(double r, int jj) if (chem_flag) jelem = element[jj]; else jelem = 0; - + double* ulist_r = ulist_r_ij[jj]; double* ulist_i = ulist_i_ij[jj]; @@ -1274,7 +1274,7 @@ void SNA::compute_duarray(double x, double y, double z, 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++) { @@ -1342,7 +1342,7 @@ double SNA::memory_usage() bytes += (double)nmax * sizeof(double); // drinnerij } if (chem_flag) bytes += (double)nmax * sizeof(int); // element - + return bytes; } /* ---------------------------------------------------------------------- */ diff --git a/src/ML-SNAP/sna.h b/src/ML-SNAP/sna.h index 126410fa00..03d77d52b4 100644 --- a/src/ML-SNAP/sna.h +++ b/src/ML-SNAP/sna.h @@ -69,7 +69,7 @@ class SNA : protected Pointers { double **dblist; // short neighbor list data - + void grow_rij(int); int nmax; // allocated size of short lists @@ -79,12 +79,12 @@ class SNA : protected Pointers { 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: @@ -125,7 +125,7 @@ class SNA : protected Pointers { double deltacg(int, int, int); void compute_ncoeff(); void compute_duarray(double, double, double, double, double, double, - double, double, int); + double, double, int); // Sets the style for the switching function // 0 = none From ec97d859b96313675461f35519ba71ac91a7b2e8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 14 Mar 2022 16:53:27 -0400 Subject: [PATCH 15/16] use tokenizer class for parsing --- src/ML-IAP/mliap_descriptor_snap.cpp | 128 ++++++++++++--------------- src/ML-IAP/mliap_model_nn.cpp | 94 ++++++++++---------- 2 files changed, 104 insertions(+), 118 deletions(-) diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index a914e523e2..dd9bf5d812 100644 --- a/src/ML-IAP/mliap_descriptor_snap.cpp +++ b/src/ML-IAP/mliap_descriptor_snap.cpp @@ -24,6 +24,7 @@ #include "mliap_data.h" #include "pair_mliap.h" #include "sna.h" +#include "tokenizer.h" #include #include @@ -402,59 +403,43 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) // 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"); + 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 || - strcmp(keywd,"rinnerelems") == 0 || - strcmp(keywd,"drinnerelems") == 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 (comm->me == 0) - utils::logmesg(lmp,"SNAP keyword {} {} ... \n", keywd, keyval); + if (comm->me == 0) utils::logmesg(lmp, "SNAP keyword {} \n", utils::trim(line)); - 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 (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 (strcmp(keywd,"rinnerelems") == 0) { - for (int ielem = 0; ielem < nelements; ielem++) { - rinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); - keyval = strtok(nullptr,"' \t\n\r\f"); - } + } else if (keywd == "rinnerelems") { + for (int ielem = 0; ielem < nelements; ielem++) + rinnerelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); rinnerflag = 1; - } else if (strcmp(keywd,"drinnerelems") == 0) { - for (int ielem = 0; ielem < nelements; ielem++) { - drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); - keyval = strtok(nullptr,"' \t\n\r\f"); - } + } else if (keywd == "drinnerelems") { + for (int ielem = 0; ielem < nelements; ielem++) + drinnerelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); rinnerflag = 1; } @@ -462,44 +447,43 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) // 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 (comm->me == 0) - utils::logmesg(lmp,"SNAP keyword {} {} \n", keywd, keyval); + auto keyval = words.next(); + if (comm->me == 0) utils::logmesg(lmp, "SNAP keyword {} {} \n", keywd, keyval); - if (strcmp(keywd,"nelems") == 0) { - 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"); + 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 = utils::numeric(FLERR,keyval,false,lmp); + } else if (keywd == "rcutfac") { + rcutfac = utils::numeric(FLERR, keyval, false, lmp); rcutfacflag = 1; - } else if (strcmp(keywd,"twojmax") == 0) { - twojmax = utils::inumeric(FLERR,keyval,false,lmp); + } else if (keywd == "twojmax") { + twojmax = utils::inumeric(FLERR, keyval, false, lmp); twojmaxflag = 1; - } else if (strcmp(keywd,"rfac0") == 0) - rfac0 = utils::numeric(FLERR,keyval,false,lmp); - else if (strcmp(keywd,"rmin0") == 0) - rmin0 = utils::numeric(FLERR,keyval,false,lmp); - else if (strcmp(keywd,"switchflag") == 0) - switchflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (strcmp(keywd,"bzeroflag") == 0) - bzeroflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (strcmp(keywd,"chemflag") == 0) - chemflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (strcmp(keywd,"bnormflag") == 0) - bnormflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (strcmp(keywd,"wselfallflag") == 0) - wselfallflag = utils::inumeric(FLERR,keyval,false,lmp); - else if (strcmp(keywd,"switchinnerflag") == 0) - switchinnerflag = utils::inumeric(FLERR,keyval,false,lmp); + } 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"); } } diff --git a/src/ML-IAP/mliap_model_nn.cpp b/src/ML-IAP/mliap_model_nn.cpp index 24fddda766..f068b07283 100644 --- a/src/ML-IAP/mliap_model_nn.cpp +++ b/src/ML-IAP/mliap_model_nn.cpp @@ -33,8 +33,8 @@ 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; @@ -73,8 +73,8 @@ void MLIAPModelNN::read_coeffs(char *coefffilename) if (comm->me == 0) { 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; @@ -111,8 +111,7 @@ 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 @@ -141,70 +140,73 @@ void MLIAPModelNN::read_coeffs(char *coefffilename) // strip comment, skip line if blank if ((ptr = strchr(line,'#'))) *ptr = '\0'; + nwords = utils::trim_and_count_words(line); if (nwords == 0) continue; + ValueTokenizer values(line, "\"' \t\n\t\f"); 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"); + 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"); 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; + 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++; + } } } } From 8e585620c30ddae46eef9a94cffe6a31ad7e0664 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 14 Mar 2022 17:14:40 -0400 Subject: [PATCH 16/16] enable and apply clang-format --- src/ML-IAP/mliap_model_nn.cpp | 330 ++++++++++++++++------------------ 1 file changed, 154 insertions(+), 176 deletions(-) diff --git a/src/ML-IAP/mliap_model_nn.cpp b/src/ML-IAP/mliap_model_nn.cpp index f068b07283..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,9 +69,9 @@ 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, + error->one(FLERR, "Cannot open MLIAPModel coeff file {}: {}", coefffilename, utils::getsyserror()); } @@ -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,13 +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; @@ -125,47 +123,53 @@ 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; ValueTokenizer values(line, "\"' \t\n\t\f"); - if (stats == 0) { // Header NET + 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"); + if (tstr.substr(0, 3) != "NET") + error->all(FLERR, "Incorrect format in MLIAPModel coefficient file"); 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 = values.next_string(); nnodes[ilayer] = values.next_int(); - 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; + 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; @@ -173,7 +177,7 @@ void MLIAPModelNN::read_coeffs(char *coefffilename) } else if (stats == 1) { scale[ielem][0][l] = values.next_double(); for (int icoeff = 1; icoeff < nwords; icoeff++) { - scale[ielem][0][l+icoeff] = values.next_double(); + scale[ielem][0][l + icoeff] = values.next_double(); } l += nwords; if (l == ndescriptors) { @@ -184,7 +188,7 @@ void MLIAPModelNN::read_coeffs(char *coefffilename) } else if (stats == 2) { scale[ielem][1][l] = values.next_double(); for (int icoeff = 1; icoeff < nwords; icoeff++) { - scale[ielem][1][l+icoeff] = values.next_double(); + scale[ielem][1][l + icoeff] = values.next_double(); } l += nwords; if (l == ndescriptors) { @@ -195,11 +199,11 @@ void MLIAPModelNN::read_coeffs(char *coefffilename) // 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] = values.next_double(); for (int icoeff = 1; icoeff < nwords; icoeff++) { - coeffelem[ielem][l+icoeff] = values.next_double(); + coeffelem[ielem][l + icoeff] = values.next_double(); } l += nwords; if (l == nparams) { @@ -216,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"); } /* ---------------------------------------------------------------------- @@ -393,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"); } /* ---------------------------------------------------------------------- @@ -410,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; }