From 11d3b13e31f48533c110ee167cf265d42c380d33 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Thu, 10 Mar 2022 18:40:59 -0700 Subject: [PATCH] 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