From 417e92bc2db164283e13fe068dbd9d7ea23d9e3a Mon Sep 17 00:00:00 2001 From: kipbarrett Date: Thu, 11 Feb 2021 08:43:04 -0600 Subject: [PATCH] Axels requested revisions --- cmake/CMakeLists.txt | 2 +- src/Makefile | 2 +- src/USER-RANN/pair_rann.cpp | 1261 +++++++++++++++++ src/USER-RANN/pair_rann.h | 174 +++ src/USER-RANN/rann_activation.h | 73 + src/USER-RANN/rann_activation_linear.h | 77 + src/USER-RANN/rann_activation_sig_i.h | 77 + src/USER-RANN/rann_fingerprint.cpp | 116 ++ src/USER-RANN/rann_fingerprint.h | 72 + src/USER-RANN/rann_fingerprint_bond.cpp | 808 +++++++++++ src/USER-RANN/rann_fingerprint_bond.h | 80 ++ .../rann_fingerprint_bondscreened.cpp | 762 ++++++++++ src/USER-RANN/rann_fingerprint_bondscreened.h | 80 ++ .../rann_fingerprint_bondscreenedspin.cpp | 814 +++++++++++ .../rann_fingerprint_bondscreenedspin.h | 80 ++ src/USER-RANN/rann_fingerprint_bondspin.cpp | 887 ++++++++++++ src/USER-RANN/rann_fingerprint_bondspin.h | 79 ++ src/USER-RANN/rann_fingerprint_radial.cpp | 239 ++++ src/USER-RANN/rann_fingerprint_radial.h | 70 + .../rann_fingerprint_radialscreened.cpp | 253 ++++ .../rann_fingerprint_radialscreened.h | 72 + .../rann_fingerprint_radialscreenedspin.cpp | 264 ++++ .../rann_fingerprint_radialscreenedspin.h | 72 + src/USER-RANN/rann_fingerprint_radialspin.cpp | 252 ++++ src/USER-RANN/rann_fingerprint_radialspin.h | 69 + src/USER-RANN/rann_list_activation.h | 2 + src/USER-RANN/rann_list_fingerprint.h | 8 + 27 files changed, 6743 insertions(+), 2 deletions(-) create mode 100644 src/USER-RANN/pair_rann.cpp create mode 100644 src/USER-RANN/pair_rann.h create mode 100644 src/USER-RANN/rann_activation.h create mode 100644 src/USER-RANN/rann_activation_linear.h create mode 100644 src/USER-RANN/rann_activation_sig_i.h create mode 100644 src/USER-RANN/rann_fingerprint.cpp create mode 100644 src/USER-RANN/rann_fingerprint.h create mode 100644 src/USER-RANN/rann_fingerprint_bond.cpp create mode 100644 src/USER-RANN/rann_fingerprint_bond.h create mode 100644 src/USER-RANN/rann_fingerprint_bondscreened.cpp create mode 100644 src/USER-RANN/rann_fingerprint_bondscreened.h create mode 100644 src/USER-RANN/rann_fingerprint_bondscreenedspin.cpp create mode 100644 src/USER-RANN/rann_fingerprint_bondscreenedspin.h create mode 100644 src/USER-RANN/rann_fingerprint_bondspin.cpp create mode 100644 src/USER-RANN/rann_fingerprint_bondspin.h create mode 100644 src/USER-RANN/rann_fingerprint_radial.cpp create mode 100644 src/USER-RANN/rann_fingerprint_radial.h create mode 100644 src/USER-RANN/rann_fingerprint_radialscreened.cpp create mode 100644 src/USER-RANN/rann_fingerprint_radialscreened.h create mode 100644 src/USER-RANN/rann_fingerprint_radialscreenedspin.cpp create mode 100644 src/USER-RANN/rann_fingerprint_radialscreenedspin.h create mode 100644 src/USER-RANN/rann_fingerprint_radialspin.cpp create mode 100644 src/USER-RANN/rann_fingerprint_radialspin.h create mode 100644 src/USER-RANN/rann_list_activation.h create mode 100644 src/USER-RANN/rann_list_fingerprint.h diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index c29fba5957..b125556c73 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -110,7 +110,7 @@ set(STANDARD_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS DIPOLE USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-MESODPD USER-CGSDK USER-COLVARS USER-DIFFRACTION USER-DPD USER-DRUDE USER-EFF USER-FEP USER-H5MD USER-LB USER-MANIFOLD USER-MEAMC USER-MESONT USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE - USER-NETCDF USER-PHONON USER-PLUMED USER-PTM USER-QTB USER-REACTION + USER-NETCDF USER-PHONON USER-PLUMED USER-PTM USER-QTB USER-RANN USER-REACTION USER-REAXC USER-SCAFACOS USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK USER-QUIP USER-QMMM USER-YAFF USER-ADIOS) set(SUFFIX_PACKAGES CORESHELL USER-OMP KOKKOS OPT USER-INTEL GPU) diff --git a/src/Makefile b/src/Makefile index 7a5e1aa728..f77e127264 100644 --- a/src/Makefile +++ b/src/Makefile @@ -56,7 +56,7 @@ PACKUSER = user-adios user-atc user-awpmd user-bocs user-cgdna user-cgsdk user-c user-intel user-lb user-manifold user-meamc user-mesodpd user-mesont \ user-mgpt user-misc user-mofff user-molfile \ user-netcdf user-omp user-phonon user-plumed user-ptm user-qmmm \ - user-qtb user-quip user-reaction user-reaxc user-scafacos user-smd user-smtbq \ + user-qtb user-quip user-rann user-reaction user-reaxc user-scafacos user-smd user-smtbq \ user-sdpd user-sph user-tally user-uef user-vtk user-yaff PACKLIB = compress gpu kim kokkos latte message mpiio mscg poems \ diff --git a/src/USER-RANN/pair_rann.cpp b/src/USER-RANN/pair_rann.cpp new file mode 100644 index 0000000000..3207b76934 --- /dev/null +++ b/src/USER-RANN/pair_rann.cpp @@ -0,0 +1,1261 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ +#define MAXLINE 1024 +#include "pair_rann.h" +#include "rann_list_activation.h" +#include "rann_list_fingerprint.h" +#include "tokenizer.h" +#include +#include +#include +#include +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "memory.h" +#include "neigh_request.h" +#include "memory.h" +#include "error.h" +#include "update.h" +#include "math_special.h" + +using namespace LAMMPS_NS; + +static const char cite_user_rann_package[] = + "USER-RANN package:\n\n" + "@Article{Nitol2021,\n" + " author = {Nitol, Mashroor S and Dickel, Doyl E and Barrett, Christopher D},\n" + " title = {Artificial neural network potential for pure zinc},\n" + " journal = {Computational Materials Science},\n" + " year = 2021,\n" + " volume = 188,\n" + " pages = {110207}\n" + "}\n\n"; + + +PairRANN::PairRANN(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + restartinfo = 0; + one_coeff = 1; + manybody_flag = 1; + allocated = 0; + nelements = -1; + elements = NULL; + mass = NULL; + + // set comm size needed by this Pair + // comm unused for now. + + comm_forward = 38; + comm_reverse = 30; + res = 10000; + cutmax = 0; + //at least one of the following will change during fingerprint definition: + doscreen = false; + allscreen = true; + + dospin = false; + fingerprint_map = new FingerprintCreatorMap(); + #define FINGERPRINT_CLASS + #define FingerprintStyle(key,Class) \ + (*fingerprint_map)[#key] = &fingerprint_creator; + + #include "rann_list_fingerprint.h" + #undef FingerprintStyle + #undef FINGERPRINT_CLASS + activation_map = new ActivationCreatorMap(); + #define ACTIVATION_CLASS + #define ActivationStyle(key,Class) \ + (*activation_map)[#key] = &activation_creator; + + #include "rann_list_activation.h" + #undef ActivationStyle + #undef ACTIVATION_CLASS +} + +PairRANN::~PairRANN() +{ + //clear memory + delete [] mass; + for (int i=0;i0) { + for (int j=0;j0) { + delete [] fingerprints[i]; + delete [] activation[i]; + } + } + delete [] fingerprints; + delete [] activation; + delete [] fingerprintcount; + delete [] fingerprintperelement; + delete [] fingerprintlength; + delete [] screening_min; + delete [] screening_max; + memory->destroy(xn); + memory->destroy(yn); + memory->destroy(zn); + memory->destroy(tn); + memory->destroy(jl); + memory->destroy(features); + memory->destroy(dfeaturesx); + memory->destroy(dfeaturesy); + memory->destroy(dfeaturesz); + memory->destroy(layer); + memory->destroy(sum); + memory->destroy(sum1); + memory->destroy(dlayerx); + memory->destroy(dlayery); + memory->destroy(dlayerz); + memory->destroy(dlayersumx); + memory->destroy(dlayersumy); + memory->destroy(dlayersumz); + if (doscreen) { + memory->destroy(Sik); + memory->destroy(Bij); + memory->destroy(dSikx); + memory->destroy(dSiky); + memory->destroy(dSikz); + memory->destroy(dSijkx); + memory->destroy(dSijky); + memory->destroy(dSijkz); + memory->destroy(dSijkxc); + memory->destroy(dSijkyc); + memory->destroy(dSijkzc); + } + if (dospin) { + memory->destroy(sx); + memory->destroy(sy); + memory->destroy(sz); + memory->destroy(dlayerx); + memory->destroy(dlayery); + memory->destroy(dlayerz); + memory->destroy(dlayersumx); + memory->destroy(dlayersumy); + memory->destroy(dlayersumz); + } +} + + + +void PairRANN::allocate(std::vector elementwords) +{ + int i,j,k,l,n; + n = atom->ntypes; + memory->create(setflag,n+1,n+1,"pair:setflag"); + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + cutmax = 0; + nmax1 = 100; + nmax2 = 20; + fmax = 0; + fnmax = 0; + nelementsp=nelements+1; + //initialize arrays + elements = new char *[nelements]; + elementsp = new char *[nelementsp];//elements + 'all' + mass = new double[nelements]; + net = new NNarchitecture[nelementsp]; + weightdefined = new bool*[nelementsp]; + biasdefined = new bool *[nelementsp]; + activation = new RANN::Activation**[nelementsp]; + fingerprints = new RANN::Fingerprint**[nelementsp]; + fingerprintlength = new int[nelementsp]; + fingerprintperelement = new int [nelementsp]; + fingerprintcount = new int[nelementsp]; + screening_min = new double [nelements*nelements*nelements]; + screening_max = new double [nelements*nelements*nelements]; + for (i=0;i 0) error->one(FLERR,"Illegal pair_style command"); +} + +void PairRANN::coeff(int narg, char **arg) +{ + int i,j; + map = new int [atom->ntypes+1]; + if (narg != 3 + atom->ntypes) error->one(FLERR,"Incorrect args for pair coefficients"); + if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) error->one(FLERR,"Incorrect args for pair coefficients"); + nelements = -1; + read_file(arg[2]); + // read args that map atom types to elements in potential file + // map[i] = which element the Ith atom type is, -1 if NULL + for (i = 3; i < narg; i++) { + if (strcmp(arg[i],"NULL") == 0) { + map[i-2] = -1; + continue; + } + for (j = 0; j < nelements; j++) { + if (strcmp(arg[i],elements[j]) == 0) break; + } + if (j < nelements) map[i-2] = j; + else error->one(FLERR,"No matching element in NN potential file"); + } + // clear setflag since coeff() called once with I,J = * * + int n = atom->ntypes; + for (i = 1; i <= n; i++) { + for (j = i; j <= n; j++) { + setflag[i][j] = 0; + } + } + // set setflag i,j for type pairs where both are mapped to elements + // set mass of atom type if i = j + int count = 0; + for (i = 1; i <= n; i++) { + for (j = i; j <= n; j++) { + if (map[i] >= 0 && map[j] >= 0) { + setflag[i][j] = 1; + if (i == j) atom->set_mass(FLERR,i,mass[map[i]]); + count++; + } + } + } + if (count == 0) error->one(FLERR,"Incorrect args for pair coefficients"); + for (i=0;iallocate(); + } + } + allocated=1; +} + +void PairRANN::read_file(char *filename) +{ + FILE *fp; + int eof = 0,i,j,k,l; + int n,nwords; + std::string line,line1; + int longline = 4096; + char linetemp[longline]; + char *ptr; + bool comment; + char str[128]; + std::vector linev,line1v; + fp = utils::open_potential(filename,lmp,nullptr); + if (fp == nullptr) {error->one(FLERR,"Cannot open RANN potential file");} + ptr=fgets(linetemp,longline,fp); + strip_comments(linetemp,fp); + line=linetemp; + while (eof == 0) { + ptr=fgets(linetemp,longline,fp); + if (ptr == NULL) { + if (check_potential()) { + error->one(FLERR,"Invalid syntax in potential file, values are inconsistent or missing"); + } + else { + eof = 1; + break; + } + } + strip_comments(linetemp,fp); + line1=linetemp; + Tokenizer values = Tokenizer(line,": ,\t_\n"); + Tokenizer values1 = Tokenizer(line1,": ,\t_\n"); + linev = values.as_vector(); + line1v = values1.as_vector(); + if (linev[0].compare("atomtypes")==0)read_atom_types(linev,line1v); + else if (linev[0].compare("mass")==0)read_mass(linev,line1v); + else if (linev[0].compare("fingerprintsperelement")==0)read_fpe(linev,line1v); + else if (linev[0].compare("fingerprints")==0)read_fingerprints(linev,line1v); + else if (linev[0].compare("fingerprintconstants")==0)read_fingerprint_constants(linev,line1v); + else if (linev[0].compare("networklayers")==0)read_network_layers(linev,line1v); + else if (linev[0].compare("layersize")==0)read_layer_size(linev,line1v); + else if (linev[0].compare("weight")==0)read_weight(linev,line1v,fp); + else if (linev[0].compare("bias")==0)read_bias(linev,line1v,fp); + else if (linev[0].compare("activationfunctions")==0)read_activation_functions(linev,line1v); + else if (linev[0].compare("calibrationparameters")==0)continue;//information on how the network was trained + else if (linev[0].compare("screening")==0)read_screening(linev,line1v); + else error->one(FLERR,"Could not understand file syntax: unknown keyword"); + ptr=fgets(linetemp,longline,fp); + strip_comments(linetemp,fp); + if (ptr == NULL) { + if (check_potential()) { + error->one(FLERR,"Invalid syntax in potential file, values are inconsistent or missing"); + } + else { + eof = 1; + break; + } + } + line=linetemp; + } +} + +void PairRANN::strip_comments(char * line,FILE* fp) { + char *ptr; + bool comment; + utils::trim_comment(line); + if (strlen(line)==0) { + comment = true; + while (comment==true) { + ptr = fgets(line,4096,fp); + if (ptr==NULL) error->one(FLERR,"Unexpected end of RANN file"); + utils::trim_comment(line); + if (strlen(line) == 0) continue; + comment = false; + } + } +} + +void PairRANN::read_atom_types(std::vector line,std::vector line1) { + int nwords = line1.size(); + if (nwords < 1) error->one(FLERR,"Incorrect syntax for atom types"); + nelements = nwords; + line1.resize(nwords+1); + line1[nwords] = "all"; + allocate(line1); +} + +void PairRANN::read_mass(std::vector line,std::vector line1) { + if (nelements == -1)error->one(FLERR,"atom types must be defined before mass in potential file."); + int nwords = 0,i; + for (i=0;ione(FLERR,"mass element not found in atom types."); +} + +void PairRANN::read_fpe(std::vector line,std::vector line1) { + int i,j; + if (nelements == -1)error->one(FLERR,"atom types must be defined before fingerprints per element in potential file."); + for (i=0;ione(FLERR,"fingerprint-per-element element not found in atom types"); +} + +void PairRANN::read_fingerprints(std::vector line,std::vector line1) { + int nwords1,nwords,i,j,k,l,m,i1; + bool found; + char str[MAXLINE]; + nwords1 = line1.size(); + nwords = line.size(); + if (nelements == -1)error->one(FLERR,"atom types must be defined before fingerprints in potential file."); + int atomtypes[nwords-1]; + for (i=1;ione(FLERR,"fingerprint element not found in atom types");} + } + i = atomtypes[0]; + k = 0; + if (fingerprintperelement[i]==-1) {error->one(FLERR,"fingerprint per element must be defined before fingerprints");} + while (kn_body_type!=nwords-1) {error->one(FLERR,"invalid fingerprint for element combination");} + k++; + fingerprints[i][i1]->init(atomtypes,strtol(line1[k++].c_str(),NULL,10)); + fingerprintcount[i]++; + } +} + + +void PairRANN::read_fingerprint_constants(std::vector line,std::vector line1) { + int i,j,k,l,m,i1; + bool found; + char str [128]; + int nwords = line.size(); + if (nelements == -1)error->one(FLERR,"atom types must be defined before fingerprints in potential file."); + int n_body_type = nwords-4; + int atomtypes[n_body_type]; + for (i=1;i<=n_body_type;i++) { + found = false; + for (j=0;jone(FLERR,"fingerprint element not found in atom types");} + } + i = atomtypes[0]; + found = false; + for (k=0;kempty) {continue;} + if (n_body_type!=fingerprints[i][k]->n_body_type) {continue;} + for (j=0;jatomtypes[j]!=atomtypes[j]) {break;} + if (j==n_body_type-1) { + if (line[nwords-3].compare(fingerprints[i][k]->style)==0 && strtol(line[nwords-2].c_str(),NULL,10)==fingerprints[i][k]->id) { + found=true; + i1 = k; + break; + } + } + } + if (found) {break;} + } + if (!found) {error->one(FLERR,"cannot define constants for unknown fingerprint");} + fingerprints[i][i1]->fullydefined=fingerprints[i][i1]->parse_values(line[nwords-1],line1); +} + +void PairRANN::read_network_layers(std::vector line,std::vector line1) { + int i,j; + if (nelements == -1)error->one(FLERR,"atom types must be defined before network layers in potential file."); + for (i=0;ione(FLERR,"invalid number of network layers"); + delete [] net[i].dimensions; + weightdefined[i] = new bool [net[i].layers]; + biasdefined[i] = new bool [net[i].layers]; + net[i].dimensions = new int [net[i].layers]; + net[i].Weights = new double * [net[i].layers-1]; + net[i].Biases = new double * [net[i].layers-1]; + net[i].activations = new int [net[i].layers-1]; + for (j=0;jone(FLERR,"network layers element not found in atom types"); +} + +void PairRANN::read_layer_size(std::vector line,std::vector line1) { + int i; + for (i=0;ione(FLERR,"networklayers for each atom type must be defined before the corresponding layer sizes."); + int j = strtol(line[2].c_str(),NULL,10); + if (j>=net[i].layers || j<0) {error->one(FLERR,"invalid layer in layer size definition");}; + net[i].dimensions[j]= strtol(line1[0].c_str(),NULL,10); + return; + } + } + error->one(FLERR,"layer size element not found in atom types"); +} + +void PairRANN::read_weight(std::vector line,std::vector line1,FILE* fp) { + int i,j,k,l,nwords; + char *ptr; + int longline = 4096; + char linetemp [longline]; + for (l=0;lone(FLERR,"networklayers must be defined before weights."); + i=strtol(line[2].c_str(),NULL,10); + if (i>=net[l].layers || i<0)error->one(FLERR,"invalid weight layer"); + if (net[l].dimensions[i]==0 || net[l].dimensions[i+1]==0) error->one(FLERR,"network layer sizes must be defined before corresponding weight"); + net[l].Weights[i] = new double [net[l].dimensions[i]*net[l].dimensions[i+1]]; + weightdefined[l][i] = true; + nwords = line1.size(); + if (nwords != net[l].dimensions[i])error->one(FLERR,"invalid weights per line"); + for (k=0;kone(FLERR,"unexpected end of potential file!"); + nwords = line1.size(); + if (nwords != net[l].dimensions[i])error->one(FLERR,"invalid weights per line"); + for (k=0;kone(FLERR,"weight element not found in atom types"); +} + +void PairRANN::read_bias(std::vector line,std::vector line1,FILE* fp) { + int i,j,l,nwords; + char *ptr; + char linetemp[MAXLINE]; + for (l=0;lone(FLERR,"networklayers must be defined before biases."); + i=strtol(line[2].c_str(),NULL,10); + if (i>=net[l].layers || i<0)error->one(FLERR,"invalid bias layer"); + if (net[l].dimensions[i]==0) error->one(FLERR,"network layer sizes must be defined before corresponding bias"); + biasdefined[l][i] = true; + net[l].Biases[i] = new double [net[l].dimensions[i+1]]; + net[l].Biases[i][0] = strtod(line1[0].c_str(),NULL); + for (j=1;jone(FLERR,"bias element not found in atom types"); +} + +void PairRANN::read_activation_functions(std::vector line,std::vector line1) { + int i,j,l,nwords; + int *ptr; + for (l=0;lone(FLERR,"networklayers must be defined before activation functions."); + i = strtol(line[2].c_str(),NULL,10); + if (i>=net[l].layers || i<0)error->one(FLERR,"invalid activation layer"); + delete activation[l][i]; + activation[l][i]=create_activation(line1[0].c_str()); + return; + } + } + error->one(FLERR,"activation function element not found in atom types"); +} + +void PairRANN::read_screening(std::vector line,std::vector line1) { + int i,j,k; + bool found; + int nwords = line.size(); + if (nelements == -1)error->one(FLERR,"atom types must be defined before fingerprints in potential file."); + if (nwords!=5)error->one(FLERR,"invalid screening command"); + int n_body_type = 3; + int atomtypes[n_body_type]; + for (i=1;i<=n_body_type;i++) { + found = false; + for (j=0;jone(FLERR,"fingerprint element not found in atom types");} + } + i = atomtypes[0]; + j = atomtypes[1]; + k = atomtypes[2]; + int index = i*nelements*nelements+j*nelements+k; + int index1 = i*nelements*nelements+k*nelements+j; + if (line[4].compare("Cmin")==0) { + screening_min[index] = strtod(line1[0].c_str(),NULL); + screening_min[index1] = screening_min[index]; + } + else if (line[4].compare("Cmax")==0) { + screening_max[index] = strtod(line1[0].c_str(),NULL); + screening_max[index1] = screening_max[index]; + } + else error->one(FLERR,"unrecognized screening keyword"); +} + +//Called after finishing reading the potential file to make sure it is complete. True is bad. +//also does the rest of the memory allocation. +bool PairRANN::check_potential() { + int i,j,k,l; + if (nelements==-1) {return true;} + for (i=0;i<=nelements;i++) { + if (inet[i].maxlayer)net[i].maxlayer = net[i].dimensions[j]; + } + if (net[i].maxlayer>fnmax) {fnmax = net[i].maxlayer;} + if (net[i].dimensions[net[i].layers-1]!=1)return true;//output layer must have single neuron (the energy) + if (net[i].dimensions[0]>fmax)fmax=net[i].dimensions[0]; + for (j=0;jempty)return true;//undefined activations + for (k=0;kfullydefined==false)return true; + fingerprints[i][j]->startingneuron = fingerprintlength[i]; + fingerprintlength[i] +=fingerprints[i][j]->get_length(); + if (fingerprints[i][j]->rc>cutmax) {cutmax = fingerprints[i][j]->rc;} + } + if (net[i].dimensions[0]!=fingerprintlength[i])return true; + } + memory->create(xn,nmax1,"pair:xn"); + memory->create(yn,nmax1,"pair:yn"); + memory->create(zn,nmax1,"pair:zn"); + memory->create(tn,nmax1,"pair:tn"); + memory->create(jl,nmax1,"pair:jl"); + memory->create(features,fmax,"pair:features"); + memory->create(dfeaturesx,fmax*nmax2,"pair:dfeaturesx"); + memory->create(dfeaturesy,fmax*nmax2,"pair:dfeaturesy"); + memory->create(dfeaturesz,fmax*nmax2,"pair:dfeaturesz"); + memory->create(layer,fnmax,"pair:layer"); + memory->create(sum,fnmax,"pair:sum"); + memory->create(sum1,fnmax,"pair:sum1"); + memory->create(dlayerx,nmax2,fnmax,"pair:dlayerx"); + memory->create(dlayery,nmax2,fnmax,"pair:dlayery"); + memory->create(dlayerz,nmax2,fnmax,"pair:dlayerz"); + memory->create(dlayersumx,nmax2,fnmax,"pair:dlayersumx"); + memory->create(dlayersumy,nmax2,fnmax,"pair:dlayersumy"); + memory->create(dlayersumz,nmax2,fnmax,"pair:dlayersumz"); + if (doscreen) { + memory->create(Sik,nmax2,"pair:Sik"); + memory->create(Bij,nmax2,"pair:Bij"); + memory->create(dSikx,nmax2,"pair:dSikx"); + memory->create(dSiky,nmax2,"pair:dSiky"); + memory->create(dSikz,nmax2,"pair:dSikz"); + memory->create(dSijkx,nmax2*nmax2,"pair:dSijkx"); + memory->create(dSijky,nmax2*nmax2,"pair:dSijky"); + memory->create(dSijkz,nmax2*nmax2,"pair:dSijkz"); + memory->create(dSijkxc,nmax2,nmax2,"pair:dSijkxc"); + memory->create(dSijkyc,nmax2,nmax2,"pair:dSijkyc"); + memory->create(dSijkzc,nmax2,nmax2,"pair:dSijkzc"); + } + if (dospin) { + memory->create(sx,fmax*nmax2,"pair:sx"); + memory->create(sy,fmax*nmax2,"pair:sy"); + memory->create(sz,fmax*nmax2,"pair:sz"); + memory->create(dlayerx,nmax2,fnmax,"pair:dsx"); + memory->create(dlayery,nmax2,fnmax,"pair:dsy"); + memory->create(dlayerz,nmax2,fnmax,"pair:dsz"); + memory->create(dlayersumx,nmax2,fnmax,"pair:dssumx"); + memory->create(dlayersumy,nmax2,fnmax,"pair:dssumy"); + memory->create(dlayersumz,nmax2,fnmax,"pair:dssumz"); + } + return false;//everything looks good +} + +void PairRANN::compute(int eflag, int vflag) +{ + //perform force/energy computation_ + if (dospin) { + if (strcmp(update->unit_style,"metal") != 0) + error->one(FLERR,"Spin pair styles require metal units"); + if (!atom->sp_flag) + error->one(FLERR,"Spin pair styles requires atom/spin style"); + } + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = vflag_atom = 0; + int ii,i,j; + int nn = 0; + sims = new Simulation[1]; + sims->inum = listfull->inum; + sims->ilist=listfull->ilist; + sims->id = listfull->ilist; + sims->type = atom->type; + sims->x = atom->x; + sims->numneigh = listfull->numneigh; + sims->firstneigh = listfull->firstneigh; + if (dospin) { + sims->s = atom->sp; + } + int itype,f,jnum,len; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = eflag_global = eflag_atom = 0; + if (eflag_global) {eng_vdwl=0;eng_coul=0;} + double energy=0; + double **force = atom->f; + double **fm = atom->fm; + double **virial = vatom; + //loop over atoms + for (ii=0;iiinum;ii++) { + i = sims->ilist[ii]; + itype = map[sims->type[i]]; + f = net[itype].dimensions[0]; + jnum = sims->numneigh[i]; + if (jnum>nmax1) { + nmax1 = jnum; + memory->grow(xn,nmax1,"pair:xn"); + memory->grow(yn,nmax1,"pair:yn"); + memory->grow(zn,nmax1,"pair:zn"); + memory->grow(tn,nmax1,"pair:tn"); + memory->grow(jl,nmax1,"pair:jl"); + } + cull_neighbor_list(&jnum,i,0); + if (jnum>nmax2) { + nmax2=jnum; + memory->grow(dfeaturesx,fmax*nmax2,"pair:dfeaturesx"); + memory->grow(dfeaturesy,fmax*nmax2,"pair:dfeaturesy"); + memory->grow(dfeaturesz,fmax*nmax2,"pair:dfeaturesz"); + memory->grow(layer,fnmax,"pair:layer"); + memory->grow(sum,fnmax,"pair:sum"); + memory->grow(sum1,fnmax,"pair:sum1"); + memory->grow(dlayerx,nmax2,fnmax,"pair:dlayerx"); + memory->grow(dlayery,nmax2,fnmax,"pair:dlayery"); + memory->grow(dlayerz,nmax2,fnmax,"pair:dlayerz"); + memory->grow(dlayersumx,nmax2,fnmax,"pair:dlayersumx"); + memory->grow(dlayersumy,nmax2,fnmax,"pair:dlayersumy"); + memory->grow(dlayersumz,nmax2,fnmax,"pair:dlayersumz"); + if (doscreen) { + memory->grow(Sik,nmax2,"pair:Sik"); + memory->grow(Bij,nmax2,"pair:Bij"); + memory->grow(dSikx,nmax2,"pair:dSikx"); + memory->grow(dSiky,nmax2,"pair:dSiky"); + memory->grow(dSikz,nmax2,"pair:dSikz"); + memory->grow(dSijkx,nmax2*nmax2,"pair:dSijkx"); + memory->grow(dSijky,nmax2*nmax2,"pair:dSijky"); + memory->grow(dSijkz,nmax2*nmax2,"pair:dSijkz"); + memory->destroy(dSijkxc); + memory->destroy(dSijkyc); + memory->destroy(dSijkzc); + memory->create(dSijkxc,nmax2,nmax2,"pair:dSijkxc"); + memory->create(dSijkyc,nmax2,nmax2,"pair:dSijkyc"); + memory->create(dSijkzc,nmax2,nmax2,"pair:dSijkzc"); + } + if (dospin) { + memory->grow(sx,fmax*nmax2,"pair:sx"); + memory->grow(sy,fmax*nmax2,"pair:sy"); + memory->grow(sz,fmax*nmax2,"pair:sz"); + memory->grow(dsx,nmax2,fnmax,"pair:dsx"); + memory->grow(dsy,nmax2,fnmax,"pair:dsy"); + memory->grow(dsz,nmax2,fnmax,"pair:dsz"); + memory->grow(dssumx,nmax2,fnmax,"pair:dssumx"); + memory->grow(dssumy,nmax2,fnmax,"pair:dssumy"); + memory->grow(dssumz,nmax2,fnmax,"pair:dssumz"); + } + } + for (j=0;jspin==false && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,ii,nn,xn,yn,zn,tn,jnum-1,jl); + else if (fingerprints[itype][j]->spin==false && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl); + else if (fingerprints[itype][j]->spin==true && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,ii,nn,xn,yn,zn,tn,jnum-1,jl); + else if (fingerprints[itype][j]->spin==true && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl); + } + itype = nelements; + //do fingerprints for type "all" + len = fingerprintperelement[itype]; + for (j=0;jspin==false && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,ii,nn,xn,yn,zn,tn,jnum-1,jl); + else if (fingerprints[itype][j]->spin==false && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl); + else if (fingerprints[itype][j]->spin==true && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,ii,nn,xn,yn,zn,tn,jnum-1,jl); + else if (fingerprints[itype][j]->spin==true && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl); + } + //run fingerprints through network + if (dospin) { + propagateforwardspin(&energy,force,fm,virial,ii,jnum); + } + else { + propagateforward(&energy,force,virial,ii,jnum); + } + } + if (vflag_fdotr) virial_fdotr_compute(); +} + +void PairRANN::cull_neighbor_list(int* jnum,int i,int sn) { + int *jlist,j,count,jj,*type,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double **x = sims[sn].x; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + type = sims[sn].type; + jlist = sims[sn].firstneigh[i]; + count = 0; + for (jj=0;jjcutmax*cutmax) { + continue; + } + xn[count]=delx; + yn[count]=dely; + zn[count]=delz; + tn[count]=jtype; + jl[count]=j; + count++; + } + jnum[0]=count+1; +} + +void PairRANN::screen_neighbor_list(int *jnum, int i,int sn) { + int jj,kk,count,count1; + count = 0; + for (jj=0;jjilist[ii]; + itype = map[sim->type[i]]; + for (int jj=0;jjcutmax*cutmax) { + Bij[kk]= false; + continue; + } + for (jj=0;jjcutmax*cutmax) { + Bij[jj] = false; + continue; + } + delx3 = delx2-delx; + dely3 = dely2-dely; + delz3 = delz2-delz; + rjk = delx3*delx3+dely3*dely3+delz3*delz3; + if (rik+rjk<=rij) {continue;}//bond angle > 90 degrees + if (rik+rij<=rjk) {continue;}//bond angle > 90 degrees + double Cmax = screening_max[itype*nelements*nelements+jtype*nelements+ktype]; + double Cmin = screening_min[itype*nelements*nelements+jtype*nelements+ktype]; + double temp1 = rij-rik+rjk; + Cn = temp1*temp1-4*rij*rjk; + //a few expanded computations provided for clarity: + //Cn = (rij-rik+rjk)*(rij-rik+rjk)-4*rij*rjk; + //Cd = (rij-rjk)*(rij-rjk)-rik*rik; + //Cijk = 1+2*(rik*rij+rik*rjk-rik*rik)/(rik*rik-(rij-rjk)*(rij-rjk)); + temp1 = rij-rjk; + Cd = temp1*temp1-rik*rik; + Cijk = Cn/Cd; + C = (Cijk-Cmin)/(Cmax-Cmin); + if (C>=1) {continue;} + else if (C<=0) { + Bij[kk]=false; + break; + } + dC = Cmax-Cmin; + dC *= dC; + dC *= dC; + temp1 = 1-C; + temp1 *= temp1; + temp1 *= temp1; + Sijk = 1-temp1; + Sijk *= Sijk; + Dij = 4*rik*(Cn+4*rjk*(rij+rik-rjk))/Cd/Cd; + Dik = -4*(rij*Cn+rjk*Cn+8*rij*rik*rjk)/Cd/Cd; + Djk = 4*rik*(Cn+4*rij*(rik-rij+rjk))/Cd/Cd; + temp1 = Cijk-Cmax; + double temp2 = temp1*temp1; + dfc = 8*temp1*temp2/(temp2*temp2-dC); + Sik[kk] *= Sijk; + dSijkx[kk*jnum+jj] = dfc*(delx*Dij-delx3*Djk); + dSikx[kk] += dfc*(delx2*Dik+delx3*Djk); + dSijky[kk*jnum+jj] = dfc*(dely*Dij-dely3*Djk); + dSiky[kk] += dfc*(dely2*Dik+dely3*Djk); + dSijkz[kk*jnum+jj] = dfc*(delz*Dij-delz3*Djk); + dSikz[kk] += dfc*(delz2*Dik+delz3*Djk); + } + } +} + + +//Called by getproperties. Propagate features and dfeatures through network. Updates force and energy +void PairRANN::propagateforward(double * energy,double **force,double **virial, int ii,int jnum) { + int i,j,k,jj,j1,itype,i1; + int *ilist,*numneigh; + ilist = listfull->ilist; + int inum = listfull->inum; + int *type = atom->type; + i1=ilist[ii]; + itype = map[type[i1]]; + NNarchitecture net1 = net[itype]; + int L = net1.layers-1; + //energy output with forces from analytical derivatives + double dsum1; + int f = net1.dimensions[0]; + for (i=0;idactivation_function(sum[j]); + sum[j] = activation[itype][i]->activation_function(sum[j]); + if (i==L-1) { + energy[j] = sum[j]; + if (eflag_atom)eatom[i1]=sum[j]; + if (eflag_global) {eng_vdwl +=sum[j];} + } + //force propagation + for (jj=0;jjilist; + int inum = listfull->inum; + int *type = atom->type; + i1=ilist[ii]; + itype = map[type[i1]]; + NNarchitecture net1 = net[itype]; + int L = net1.layers-1; + //energy output with forces from analytical derivatives + double dsum1; + int f = net1.dimensions[0]; + for (i=0;idactivation_function(sum[j]); + sum[j] = activation[itype][i]->activation_function(sum[j]); + if (i==L-1) { + energy[j] = sum[j]; + if (eflag_atom)eatom[i1]=sum[j]; + if (eflag_global) {eng_vdwl +=sum[j];} + } + //force propagation + for (jj=0;jjrequest(this,instance_me); + neighbor->requests[irequest_full]->id = 1; + neighbor->requests[irequest_full]->half = 0; + neighbor->requests[irequest_full]->full = 1; +} + + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairRANN::init_one(int i, int j) +{ + return cutmax; +} + +void PairRANN::errorf(const char *file, int line, const char * message) { + error->one(file,line,message); +} + +int PairRANN::factorial(int n) { + return round(MathSpecial::factorial(n)); +} + +template +RANN::Fingerprint *PairRANN::fingerprint_creator(PairRANN* pair) +{ + return new T(pair); +} + +RANN::Fingerprint *PairRANN::create_fingerprint(const char *style) +{ + if (fingerprint_map->find(style) != fingerprint_map->end()) { + FingerprintCreator fingerprint_creator = (*fingerprint_map)[style]; + return fingerprint_creator(this); + } + char str[128]; + sprintf(str,"Unknown fingerprint style %s",style); + error->one(FLERR,str); + return NULL; +} + +template +RANN::Activation *PairRANN::activation_creator(PairRANN* pair) +{ + return new T(pair); +} + +RANN::Activation *PairRANN::create_activation(const char *style) +{ + if (activation_map->find(style) != activation_map->end()) { + ActivationCreator activation_creator = (*activation_map)[style]; + return activation_creator(this); + } + char str[128]; + sprintf(str,"Unknown activation style %s",style); + error->one(FLERR,str); + return NULL; +} + diff --git a/src/USER-RANN/pair_rann.h b/src/USER-RANN/pair_rann.h new file mode 100644 index 0000000000..9f1dbe9212 --- /dev/null +++ b/src/USER-RANN/pair_rann.h @@ -0,0 +1,174 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ + +#ifdef PAIR_CLASS + +PairStyle(rann,PairRANN) + +#else + +#ifndef LMP_PAIR_RANN +#define LMP_PAIR_RANN + + +#include "neighbor.h" +#include "neigh_list.h" +#include "pair.h" +#include + + +namespace LAMMPS_NS { + + namespace RANN { + //forward declarations + class Activation; + class Fingerprint; + } + + class PairRANN : public Pair { + public: + //inherited functions + PairRANN(class LAMMPS *); + ~PairRANN(); + void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + void init_style(); + double init_one(int, int); + void init_list(int , NeighList *); + void errorf(const char*,int,const char*); + int factorial(int); + //black magic for modular fingerprints and activations + class RANN::Activation ***activation; + class RANN::Fingerprint ***fingerprints; + typedef RANN::Fingerprint *(*FingerprintCreator)(PairRANN *); + typedef RANN::Activation *(*ActivationCreator)(PairRANN *); + typedef std::map FingerprintCreatorMap; + typedef std::map ActivationCreatorMap; + FingerprintCreatorMap *fingerprint_map; + ActivationCreatorMap *activation_map; + RANN::Fingerprint * create_fingerprint(const char *); + RANN::Activation * create_activation(const char *); + + //global variables + int nelements; // # of elements (distinct from LAMMPS atom types since multiple atom types can be mapped to one element) + int nelementsp; // nelements+1 + char **elements; // names of elements + char **elementsp; // names of elements with "all" appended as the last "element" + double *mass; // mass of each element + double cutmax; // max radial distance for neighbor lists + int *map; // mapping from atom types to elements + int *fingerprintcount; // static variable used in initialization + int *fingerprintlength; // # of input neurons defined by fingerprints of each element. + int *fingerprintperelement; // # of fingerprints for each element + bool doscreen;//screening is calculated if any defined fingerprint uses it + bool allscreen;//all fingerprints use screening so screened neighbors can be completely ignored + bool dospin; + int res;//Resolution of function tables for cubic interpolation. + int memguess; + double *screening_min; + double *screening_max; + bool **weightdefined; + bool **biasdefined; + int nmax1; + int nmax2; + int fmax; + int fnmax; + //memory actively written to during each compute: + double *xn,*yn,*zn,*Sik,*dSikx,*dSiky,*dSikz,*dSijkx,*dSijky,*dSijkz,*sx,*sy,*sz,**dSijkxc,**dSijkyc,**dSijkzc,*dfeaturesx,*dfeaturesy,*dfeaturesz,*features; + double *layer,*sum,*sum1,**dlayerx,**dlayery,**dlayerz,**dlayersumx,**dlayersumy,**dlayersumz; + double **dsx,**dsy,**dsz,**dssumx,**dssumy,**dssumz; + int *tn,*jl; + bool *Bij; + + struct Simulation{ + int *id; + bool forces; + bool spins; + double **x; + double **f; + double **s; + double box[3][3]; + double origin[3]; + double **features; + double **dfx; + double **dfy; + double **dfz; + double **dsx; + double **dsy; + double **dsz; + int *ilist,*numneigh,**firstneigh,*type,inum,gnum; + }; + + struct NNarchitecture{ + int layers; + int *dimensions;//vector of length layers with entries for neurons per layer + double **Weights; + double **Biases; + int *activations;//unused + int maxlayer;//longest layer (for memory allocation) + }; + + Simulation *sims; + NNarchitecture *net;//array of networks, 1 for each element. + + private: + template static RANN::Fingerprint *fingerprint_creator(PairRANN *); + template static RANN::Activation *activation_creator(PairRANN *); + //new functions + void allocate(std::vector);//called after reading element list, but before reading the rest of the potential + void read_file(char *);//read potential file + void strip_comments(char *,FILE *fp); + void read_atom_types(std::vector,std::vector); + void read_mass(std::vector,std::vector); + void read_fpe(std::vector,std::vector);//fingerprints per element. Count total fingerprints defined for each 1st element in element combinations + void read_fingerprints(std::vector,std::vector); + void read_fingerprint_constants(std::vector,std::vector); + void read_network_layers(std::vector,std::vector);//include input and output layer (hidden layers + 2) + void read_layer_size(std::vector,std::vector); + void read_weight(std::vector,std::vector,FILE*);//weights should be formatted as properly shaped matrices + void read_bias(std::vector,std::vector,FILE*);//biases should be formatted as properly shaped vectors + void read_activation_functions(std::vector,std::vector); + void read_screening(std::vector,std::vector); + bool check_potential();//after finishing reading potential file + void propagateforward(double *,double **,double **,int,int);//called by compute to get force and energy + void propagateforwardspin(double *,double **,double **,double**,int,int);//called by compute to get force and energy + void screen(int,int,int); + void cull_neighbor_list(int *,int,int); + void screen_neighbor_list(int *,int,int); + }; + +} + +#endif +#endif + + + diff --git a/src/USER-RANN/rann_activation.h b/src/USER-RANN/rann_activation.h new file mode 100644 index 0000000000..d8d2940228 --- /dev/null +++ b/src/USER-RANN/rann_activation.h @@ -0,0 +1,73 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ + +#ifndef LMP_RANN_ACTIVATION_H +#define LMP_RANN_ACTIVATION_H + + +namespace LAMMPS_NS { + namespace RANN { + class Activation { + public: + Activation(class PairRANN *); + virtual ~Activation(); + virtual double activation_function(double); + virtual double dactivation_function(double); + virtual double ddactivation_function(double); + bool empty; + const char *style; + }; + + Activation::Activation(PairRANN *_pair) { + empty = true; + style = "empty"; + } + + Activation::~Activation() { + + } + + //default is linear activation (no change). + double Activation::activation_function(double A) { + return A; + } + + double Activation::dactivation_function(double A) { + return 1.0; + } + + double Activation::ddactivation_function(double A) { + return 0.0; + } + } +} + + +#endif /* RANN_ACTIVATION_H_ */ diff --git a/src/USER-RANN/rann_activation_linear.h b/src/USER-RANN/rann_activation_linear.h new file mode 100644 index 0000000000..812c570d3a --- /dev/null +++ b/src/USER-RANN/rann_activation_linear.h @@ -0,0 +1,77 @@ + +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ +#ifdef ACTIVATION_CLASS + +ActivationStyle(linear,Activation_linear) + +#else + +#ifndef LMP_RANN_ACTIVATION_LINEAR_H +#define LMP_RANN_ACTIVATION_LINEAR_H + +#include "rann_activation.h" + +namespace LAMMPS_NS { + namespace RANN { + + class Activation_linear : public Activation { + public: + Activation_linear(class PairRANN *); + double activation_function(double); + double dactivation_function(double); + double ddactivation_function(double); + }; + + Activation_linear::Activation_linear(PairRANN *_pair) : Activation(_pair) { + empty = false; + style = "linear"; + } + + double Activation_linear::activation_function(double A) + { + return A; + } + + double Activation_linear::dactivation_function(double A) + { + return 1.0; + } + + double Activation_linear::ddactivation_function(double) { + return 0.0; + } + + + } +} + +#endif +#endif /* ACTIVATION_LINEAR_H_ */ diff --git a/src/USER-RANN/rann_activation_sig_i.h b/src/USER-RANN/rann_activation_sig_i.h new file mode 100644 index 0000000000..9ca9c55d3d --- /dev/null +++ b/src/USER-RANN/rann_activation_sig_i.h @@ -0,0 +1,77 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ + +#ifdef ACTIVATION_CLASS + +ActivationStyle(sigI,Activation_sigI) + +#else + +#ifndef LMP_RANN_ACTIVATION_SIGI_H +#define LMP_RANN_ACTIVATION_SIGI_H + +#include "rann_activation.h" + +namespace LAMMPS_NS { + namespace RANN { + + class Activation_sigI : public Activation { + public: + Activation_sigI(class PairRANN *); + double activation_function(double); + double dactivation_function(double); + double ddactivation_function(double); + }; + + Activation_sigI::Activation_sigI(PairRANN *_pair) : Activation(_pair) { + empty = false; + style = "sigI"; + } + + double Activation_sigI::activation_function(double in) { + if (in>34)return in; + return 0.1*in + 0.9*log(exp(in) + 1); + } + + double Activation_sigI::dactivation_function(double in) { + if (in>34)return 1; + return 0.1 + 0.9/(exp(in)+1)*exp(in); + } + + double Activation_sigI::ddactivation_function(double in) { + if (in>34)return 0; + return 0.9*exp(in)/(exp(in)+1)/(exp(in)+1);; + } + + } +} + +#endif +#endif /* ACTIVATION_SIGI_H_ */ diff --git a/src/USER-RANN/rann_fingerprint.cpp b/src/USER-RANN/rann_fingerprint.cpp new file mode 100644 index 0000000000..f0bf62ed62 --- /dev/null +++ b/src/USER-RANN/rann_fingerprint.cpp @@ -0,0 +1,116 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ + +#include "rann_fingerprint.h" + +using namespace LAMMPS_NS::RANN; + +Fingerprint::Fingerprint(PairRANN *_pair) +{ + spin = false; + screen = false; + empty = true; + fullydefined = false; + n_body_type = 0; + style = "empty"; + pair = _pair; +} + +Fingerprint::~Fingerprint() { + +} + +bool Fingerprint::parse_values(std::string line,std::vector line1) { + return false; +} + +void Fingerprint::init(int *i,int id) { + +} + +void Fingerprint::allocate() { + +} + +void Fingerprint::compute_fingerprint(double *features,double *dfeaturesx,double *dfeaturesy,double * dfeaturesz, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + +} + +void Fingerprint::compute_fingerprint(double *features,double *dfeaturesx,double *dfeaturesy,double * dfeaturesz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + +} + +void Fingerprint::compute_fingerprint(double *features,double *dfeaturesx,double *dfeaturesy,double * dfeaturesz,double *sx, double *sy, double *sz, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + +} + +void Fingerprint::compute_fingerprint(double *features,double *dfeaturesx,double *dfeaturesy,double * dfeaturesz,double *sx, double *sy, double *sz, double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + +} + +void Fingerprint::write_values(FILE *fid) { + +} + +int Fingerprint::get_length() { + return 0; +} + +//Smooth cutoff, goes from 1 to zero over the interval rc-dr to rc. Same as MEAM uses. Used by generateradialtable and generatexpcuttable. +double Fingerprint::cutofffunction(double r,double rc, double dr) { + double out; + if (r < (rc -dr))out=1; + else if (r>rc)out=0; + else { + out = 1-(rc-r)/dr; + out *= out; + out *= out; + out = 1-out; + out *= out; + } + return out; +} + +void Fingerprint::generate_rinvssqrttable() { + int buf = 5; + int m; + double r1; + double cutmax = pair->cutmax; + int res = pair->res; + rinvsqrttable = new double[res+buf]; + for (m=0;m<(res+buf);m++) { + r1 = cutmax*cutmax*(double)(m)/(double)(res); + rinvsqrttable[m] = 1/sqrt(r1); + } +} + + + + diff --git a/src/USER-RANN/rann_fingerprint.h b/src/USER-RANN/rann_fingerprint.h new file mode 100644 index 0000000000..8b777ce398 --- /dev/null +++ b/src/USER-RANN/rann_fingerprint.h @@ -0,0 +1,72 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + +----------------*/ + +#ifndef LMP_RANN_FINGERPRINT_H +#define LMP_RANN_FINGERPRINT_H + +#include "pair_rann.h" + +namespace LAMMPS_NS { + namespace RANN { + class Fingerprint { + public: + Fingerprint(PairRANN *); + virtual ~Fingerprint(); + virtual bool parse_values(std::string,std::vector); + virtual void write_values(FILE *); + virtual void init(int*,int); + virtual void allocate(); + void init_screen(int); + virtual void compute_fingerprint(double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*);//no screen,no spin + virtual void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);//screen + virtual void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*);//spin + virtual void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);//spin,screen + virtual int get_length(); + virtual double cutofffunction(double,double, double); + virtual void generate_rinvssqrttable(); + bool spin; + bool screen; + int n_body_type;//i-j vs. i-j-k vs. i-j-k-l, etc. + bool empty; + bool fullydefined; + int startingneuron; + int id;//based on ordering of fingerprints listed for i-j in potential file + const char *style; + int* atomtypes; + double *rinvsqrttable; + double rc; + PairRANN *pair; + }; + } +} + + +#endif /* RANN_FINGERPRINT_H_ */ diff --git a/src/USER-RANN/rann_fingerprint_bond.cpp b/src/USER-RANN/rann_fingerprint_bond.cpp new file mode 100644 index 0000000000..672cea718d --- /dev/null +++ b/src/USER-RANN/rann_fingerprint_bond.cpp @@ -0,0 +1,808 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ + +#include "rann_fingerprint_bond.h" + + +using namespace LAMMPS_NS::RANN; + +Fingerprint_bond::Fingerprint_bond(PairRANN *_pair) : Fingerprint(_pair) +{ + n_body_type = 3; + dr = 0; + re = 0; + rc = 0; + alpha_k = new double [1]; + alpha_k[0] = -1; + kmax = 0; + mlength = 0; + id = -1; + style = "bond"; + atomtypes = new int [n_body_type]; + empty = true; + _pair->allscreen = false; +} + +Fingerprint_bond::~Fingerprint_bond() { + delete [] alpha_k; + delete [] atomtypes; + delete [] expcuttable; + delete [] dfctable; + for (int i=0;i<(mlength*(mlength+1))>>1;i++) { + delete [] coeff[i]; + delete [] coeffx[i]; + delete [] coeffy[i]; + delete [] coeffz[i]; + delete [] Mf[i]; + } + delete [] coeff; + delete [] coeffx; + delete [] coeffy; + delete [] coeffz; + delete [] Mf; + delete [] rinvsqrttable; +} + +bool Fingerprint_bond::parse_values(std::string constant,std::vector line1) { + int nwords,l; + nwords=line1.size(); + if (constant.compare("re")==0) { + re = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("rc")==0) { + rc = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("alphak")==0) { + delete [] alpha_k; + alpha_k = new double [nwords]; + for (l=0;lerrorf(FLERR,"Undefined value for bond power"); + if (re!=0.0 && rc!=0.0 && alpha_k[0]!=-1 && dr!=0.0 && mlength!=0 && kmax!=0)return true; + return false; +} + +void Fingerprint_bond::write_values(FILE *fid) { + int i; + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:re:\n",style,id); + fprintf(fid,"%f\n",re); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:rc:\n",style,id); + fprintf(fid,"%f\n",rc); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:alphak:\n",style,id); + for (i=0;ielementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:dr:\n",style,id); + fprintf(fid,"%f\n",dr); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:k:\n",style,id); + fprintf(fid,"%d\n",kmax); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:m:\n",style,id); + fprintf(fid,"%d\n",mlength); +} + +void Fingerprint_bond::init(int *i,int _id) { + for (int j=0;jres; + double cutmax = pair->cutmax; + expcuttable = new double [(res+buf)*(kmax)]; + dfctable = new double [res+buf]; + for (m=0;m<(res+buf);m++) { + r1 = cutmax*cutmax*(double)(m)/(double)(res); + for (n=0;n<(kmax);n++) { + expcuttable[n+m*(kmax)] = exp(-alpha_k[n]/re*sqrt(r1))*cutofffunction(sqrt(r1),rc,dr); + } + if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) { + dfctable[m]=0; + } + else{ + dfctable[m]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4)); + } + } +} + +//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop. +void Fingerprint_bond::generate_coefficients() { //calculates multinomial coefficient for each term + int p,mb,mc; + int n,p1,i1; + mb = mlength; + mc=(mb*(mb+1))>>1; + coeff = new int *[mc]; + coeffx = new int *[mc]; + coeffy = new int *[mc]; + coeffz = new int *[mc]; + for (p=0;pfactorial(p)/pair->factorial(coeffx[p1][p])/pair->factorial(coeffy[p1][p])/pair->factorial(coeffz[p1][p]); + } + } +} + + +//Called by getproperties. Gets 3-body features and dfeatures +void Fingerprint_bond::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + int i; + int *ilist,*numneigh; + PairRANN::Simulation *sim = &pair->sims[sid]; + ilist = sim->ilist; + numneigh = sim->numneigh; + i = ilist[ii]; + //select the more efficient algorithm for this particular potential and environment. + if (jnum*2>(mlength+1)*mlength*20) { + do3bodyfeatureset_singleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,ii,sid,xn,yn,zn,tn,jnum,jl); + } + else{ + do3bodyfeatureset_doubleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,ii,sid,xn,yn,zn,tn,jnum,jl); + + } +} + +//Called by do3bodyfeatureset. Algorithm for high neighbor numbers and small series of bond angle powers +void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + int i,j,jj,itype,jtype,kk,m,n,mcount,a,a1,a2,ai; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *ilist,*jlist,*numneigh,**firstneigh; + int count=0; + PairRANN::Simulation *sim = &pair->sims[sid]; + int *type = sim->type; + double cutmax = pair->cutmax; + int res = pair->res; + double cutinv2 = 1/cutmax/cutmax; + ilist = sim->ilist; + int nelements=pair->nelements; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double expr[jnum][kmax+12]; + int p = kmax; + int countmb=((mlength)*(mlength+1))>>1; + // calculate interpolation expr, rinvs and dfc, for each neighbor + for (jj = 0; jj < jnum; jj++) { + jtype = tn[jj]; + if (atomtypes[1] != nelements && atomtypes[1] != jtype && atomtypes[2] != nelements && atomtypes[2] != jtype) { + expr[jj][0]=0; + continue; + } + delx = xn[jj]; + dely = yn[jj]; + delz = zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq>rc*rc) { + expr[jj][0]=0; + continue; + } + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + r1 = r1-trunc(r1); + double *p0 = &expcuttable[(m1-1)*kmax]; + double *p1 = &expcuttable[m1*kmax]; + double *p2 = &expcuttable[(m1+1)*kmax]; + double *p3 = &expcuttable[(m1+2)*kmax]; + for (kk=0;kksims[sid]; + double **x = sim->x; + int *type = sim->type; + int nelements = pair->nelements; + int res = pair->res; + double cutmax = pair->cutmax; + double cutinv2 = 1/cutmax/cutmax; + ilist = sim->ilist; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double expr[jnum][kmax]; + double y[jnum][3]; + double ri[jnum]; + double dfc[jnum]; + int kb = kmax; + int mb = mlength; + double c41[kmax]; + double c51[kmax]; + double c61[kmax]; + double ct[kmax]; + for (jj = 0; jj < jnum; jj++) { + jtype = tn[jj]; + if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype) { + expr[jj][0]=0; + continue; + } + delx = xn[jj]; + dely = yn[jj]; + delz = zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq>rc*rc) { + expr[jj][0]=0; + continue; + } + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + if (!(m1>=1 && m1 <= res))pair->errorf(FLERR,"Neighbor list is invalid.");//usually results from nan somewhere. + r1 = r1-trunc(r1); + double *p0 = &expcuttable[(m1-1)*kmax]; + double *p1 = &expcuttable[m1*kmax]; + double *p2 = &expcuttable[(m1+1)*kmax]; + double *p3 = &expcuttable[(m1+2)*kmax]; + for (kk=0;kk); + void write_values(FILE *); + void init(int*,int); + void allocate(); + void compute_fingerprint(double*,double*,double*,double*,int,int,double *,double *,double *,int *,int,int *); + void do3bodyfeatureset_doubleneighborloop(double *,double *,double *,double *,int,int,double *,double *,double *,int *,int,int *); + void do3bodyfeatureset_singleneighborloop(double *,double *,double *,double *,int,int,double *,double *,double *,int *,int,int *); + void generate_exp_cut_table(); + void generate_coefficients(); + int get_length(); + + double *expcuttable; + double *dfctable; + double dr; + double *alpha_k; + double re; + int **coeff; + int **coeffx; + int **coeffy; + int **coeffz; + int kmax; + int mlength; + int **Mf; + + }; + + } +} + + +#endif +#endif /* FINGERPRINT_BOND_H_ */ diff --git a/src/USER-RANN/rann_fingerprint_bondscreened.cpp b/src/USER-RANN/rann_fingerprint_bondscreened.cpp new file mode 100644 index 0000000000..5a9a9b7a4c --- /dev/null +++ b/src/USER-RANN/rann_fingerprint_bondscreened.cpp @@ -0,0 +1,762 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ + +#include "rann_fingerprint_bondscreened.h" + +using namespace LAMMPS_NS::RANN; + +Fingerprint_bondscreened::Fingerprint_bondscreened(PairRANN *_pair) : Fingerprint(_pair) +{ + n_body_type = 3; + dr = 0; + re = 0; + rc = 0; + alpha_k = new double [1]; + alpha_k[0] = -1; + kmax = 0; + mlength = 0; + id = -1; + style = "bondscreened"; + atomtypes = new int [n_body_type]; + empty = true; + _pair->doscreen = true; + screen = true; +} + +Fingerprint_bondscreened::~Fingerprint_bondscreened() { + delete [] alpha_k; + delete [] atomtypes; + delete [] expcuttable; + delete [] dfctable; + for (int i=0;i<(mlength*(mlength+1))>>1;i++) { + delete [] coeff[i]; + delete [] coeffx[i]; + delete [] coeffy[i]; + delete [] coeffz[i]; + delete [] Mf[i]; + } + delete [] coeff; + delete [] coeffx; + delete [] coeffy; + delete [] coeffz; + delete [] Mf; + delete [] rinvsqrttable; +} + +bool Fingerprint_bondscreened::parse_values(std::string constant,std::vector line1) { + int nwords,l; + nwords=line1.size(); + if (constant.compare("re")==0) { + re = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("rc")==0) { + rc = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("alphak")==0) { + delete [] alpha_k; + alpha_k = new double [nwords]; + for (l=0;lerrorf(FLERR,"Undefined value for bond power"); + if (re!=0.0 && rc!=0.0 && alpha_k[0]!=-1 && dr!=0.0 && mlength!=0 && kmax!=0)return true; + return false; +} + +void Fingerprint_bondscreened::write_values(FILE *fid) { + int i; + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:re:\n",style,id); + fprintf(fid,"%f\n",re); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:rc:\n",style,id); + fprintf(fid,"%f\n",rc); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:alphak:\n",style,id); + for (i=0;ielementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:dr:\n",style,id); + fprintf(fid,"%f\n",dr); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:k:\n",style,id); + fprintf(fid,"%d\n",kmax); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:m:\n",style,id); + fprintf(fid,"%d\n",mlength); +} + +void Fingerprint_bondscreened::init(int *i,int _id) { + for (int j=0;jres; + double cutmax = pair->cutmax; + expcuttable = new double [(res+buf)*(kmax)]; + dfctable = new double [res+buf]; + for (m=0;m<(res+buf);m++) { + r1 = cutmax*cutmax*(double)(m)/(double)(res); + for (n=0;n<(kmax);n++) { + expcuttable[n+m*(kmax)] = exp(-alpha_k[n]/re*sqrt(r1))*cutofffunction(sqrt(r1),rc,dr); + } + if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) { + dfctable[m]=0; + } + else{ + dfctable[m]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4)); + } + } +} + + +//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop. +void Fingerprint_bondscreened::generate_coefficients() { //calculates multinomial coefficient for each term + int p,mb,mc; + int n,p1,i1; + mb = mlength; + mc=(mb*(mb+1))>>1; + coeff = new int *[mc]; + coeffx = new int *[mc]; + coeffy = new int *[mc]; + coeffz = new int *[mc]; + for (p=0;pfactorial(p)/pair->factorial(coeffx[p1][p])/pair->factorial(coeffy[p1][p])/pair->factorial(coeffz[p1][p]); + } + } +} + + +//Called by getproperties. Gets 3-body features and dfeatures +void Fingerprint_bondscreened::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + int i; + //select the more efficient algorithm for this particular potential and environment. + if (jnum*2>(mlength+1)*mlength*20) { + do3bodyfeatureset_singleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,sid,xn,yn,zn,tn,jnum,jl); + } + else{ + do3bodyfeatureset_doubleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,sid,xn,yn,zn,tn,jnum,jl); + + } +} + +//Called by do3bodyfeatureset. Algorithm for high neighbor numbers and small series of bond angle powers +void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + int i,j,jj,itype,jtype,kk,m,n,mcount,a,a1,a2,ai; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *ilist,*jlist,*numneigh,**firstneigh; + int count=0; + PairRANN::Simulation *sim = &pair->sims[sid]; + int *type = sim->type; + double cutmax = pair->cutmax; + int res = pair->res; + double cutinv2 = 1/cutmax/cutmax; + ilist = sim->ilist; + int nelements=pair->nelements; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double expr[jnum][kmax+12]; + int p = kmax; + int countmb=((mlength)*(mlength+1))>>1; + // calculate interpolation expr, rinvs and dfc, for each neighbor + for (jj = 0; jj < jnum; jj++) { + if (Bij[jj]==false) {continue;} + jtype = tn[jj]; + if (atomtypes[1] != nelements && atomtypes[1] != jtype && atomtypes[2] != nelements && atomtypes[2] != jtype) { + expr[jj][0]=0; + continue; + } + delx=xn[jj]; + dely=yn[jj]; + delz=zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq>rc*rc) { + expr[jj][0]=0; + continue; + } + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + r1 = r1-trunc(r1); + double *p0 = &expcuttable[(m1-1)*kmax]; + double *p1 = &expcuttable[m1*kmax]; + double *p2 = &expcuttable[(m1+1)*kmax]; + double *p3 = &expcuttable[(m1+2)*kmax]; + for (kk=0;kksims[sid]; + int *type = sim->type; + int nelements = pair->nelements; + int res = pair->res; + double cutmax = pair->cutmax; + double cutinv2 = 1/cutmax/cutmax; + ilist = sim->ilist; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double expr[jnum][kmax]; + double y[jnum][3]; + double ri[jnum]; + double dfc[jnum]; + int kb = kmax; + int mb = mlength; + double Bijk[kb][mb]; + double c41[kmax]; + double c51[kmax]; + double c61[kmax]; + double ct[kmax]; + for (jj = 0; jj < jnum; jj++) { + if (Bij[jj]==false) {continue;} + jtype = tn[jj]; + if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype) { + expr[jj][0]=0; + continue; + } + delx = xn[jj]; + dely = yn[jj]; + delz = zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq>rc*rc) { + expr[jj][0]=0; + continue; + } + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + if (!(m1>=1 && m1 <= res))pair->errorf(FLERR,"Neighbor list is invalid.");//usually results from nan somewhere. + r1 = r1-trunc(r1); + double *p0 = &expcuttable[(m1-1)*kmax]; + double *p1 = &expcuttable[m1*kmax]; + double *p2 = &expcuttable[(m1+1)*kmax]; + double *p3 = &expcuttable[(m1+2)*kmax]; + for (kk=0;kk); + void write_values(FILE *); + void init(int*,int); + void allocate(); + void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*); + void do3bodyfeatureset_doubleneighborloop(double *,double *,double *,double *,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*); + void do3bodyfeatureset_singleneighborloop(double *,double *,double *,double *,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*); + void generate_exp_cut_table(); + void generate_coefficients(); + int get_length(); + + double *expcuttable; + double *dfctable; + double dr; + double *alpha_k; + double re; + int **coeff; + int **coeffx; + int **coeffy; + int **coeffz; + int kmax; + int mlength; + int **Mf; + + }; + } + +} + + +#endif +#endif /* FINGERPRINT_BOND_H_ */ diff --git a/src/USER-RANN/rann_fingerprint_bondscreenedspin.cpp b/src/USER-RANN/rann_fingerprint_bondscreenedspin.cpp new file mode 100644 index 0000000000..f41dde883d --- /dev/null +++ b/src/USER-RANN/rann_fingerprint_bondscreenedspin.cpp @@ -0,0 +1,814 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ +#include "rann_fingerprint_bondscreenedspin.h" + + +using namespace LAMMPS_NS::RANN; + +Fingerprint_bondscreenedspin::Fingerprint_bondscreenedspin(PairRANN *_pair) : Fingerprint(_pair) +{ + n_body_type = 3; + dr = 0; + re = 0; + rc = 0; + alpha_k = new double [1]; + alpha_k[0] = -1; + kmax = 0; + mlength = 0; + id = -1; + style = "bondscreenedspin"; + atomtypes = new int [n_body_type]; + empty = true; + _pair->doscreen = true; + screen = true; + _pair->dospin = true; + spin = true; +} + +Fingerprint_bondscreenedspin::~Fingerprint_bondscreenedspin() { + delete [] alpha_k; + delete [] atomtypes; + delete [] expcuttable; + delete [] dfctable; + for (int i=0;i<(mlength*(mlength+1))>>1;i++) { + delete [] coeff[i]; + delete [] coeffx[i]; + delete [] coeffy[i]; + delete [] coeffz[i]; + delete [] Mf[i]; + } + delete [] coeff; + delete [] coeffx; + delete [] coeffy; + delete [] coeffz; + delete [] Mf; + delete [] rinvsqrttable; +} + +bool Fingerprint_bondscreenedspin::parse_values(std::string constant,std::vector line1) { + int nwords,l; + nwords=line1.size(); + if (constant.compare("re")==0) { + re = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("rc")==0) { + rc = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("alphak")==0) { + delete [] alpha_k; + alpha_k = new double [nwords]; + for (l=0;lerrorf(FLERR,"Undefined value for bond power"); + if (re!=0.0 && rc!=0.0 && alpha_k[0]!=-1 && dr!=0.0 && mlength!=0 && kmax!=0)return true; + return false; +} + +void Fingerprint_bondscreenedspin::write_values(FILE *fid) { + int i; + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:re:\n",style,id); + fprintf(fid,"%f\n",re); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:rc:\n",style,id); + fprintf(fid,"%f\n",rc); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:alphak:\n",style,id); + for (i=0;ielementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:dr:\n",style,id); + fprintf(fid,"%f\n",dr); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:k:\n",style,id); + fprintf(fid,"%d\n",kmax); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:m:\n",style,id); + fprintf(fid,"%d\n",mlength); +} + +void Fingerprint_bondscreenedspin::init(int *i,int _id) { + for (int j=0;jres; + double cutmax = pair->cutmax; + expcuttable = new double [(res+buf)*(kmax)]; + dfctable = new double [res+buf]; + for (m=0;m<(res+buf);m++) { + r1 = cutmax*cutmax*(double)(m)/(double)(res); + for (n=0;n<(kmax);n++) { + expcuttable[n+m*(kmax)] = exp(-alpha_k[n]/re*sqrt(r1))*cutofffunction(sqrt(r1),rc,dr); + } + if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) { + dfctable[m]=0; + } + else{ + dfctable[m]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4)); + } + } +} + + +//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop. +void Fingerprint_bondscreenedspin::generate_coefficients() { //calculates multinomial coefficient for each term + int p,mb,mc; + int n,p1,i1; + mb = mlength; + mc=(mb*(mb+1))>>1; + coeff = new int *[mc]; + coeffx = new int *[mc]; + coeffy = new int *[mc]; + coeffz = new int *[mc]; + for (p=0;pfactorial(p)/pair->factorial(coeffx[p1][p])/pair->factorial(coeffy[p1][p])/pair->factorial(coeffz[p1][p]); + } + } +} + + +//Called by getproperties. Gets 3-body features and dfeatures +void Fingerprint_bondscreenedspin::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz, double *dspinx, double *dspiny, double *dspinz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + int i; + //select the more efficient algorithm for this particular potential and environment. + if (jnum*2>(mlength+1)*mlength*20) { + do3bodyfeatureset_singleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,dspinx,dspiny,dspinz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,sid,xn,yn,zn,tn,jnum,jl); + } + else{ + do3bodyfeatureset_doubleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,dspinx,dspiny,dspinz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,sid,xn,yn,zn,tn,jnum,jl); + + } +} + +//Called by do3bodyfeatureset. Algorithm for high neighbor numbers and small series of bond angle powers +void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz, double *dspinx, double *dspiny, double *dspinz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + int i,j,jj,itype,jtype,kk,m,n,mcount,a,a1,a2,ai; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *ilist,*jlist,*numneigh,**firstneigh; + int count=0; + PairRANN::Simulation *sim = &pair->sims[sid]; + int *type = sim->type; + double cutmax = pair->cutmax; + int res = pair->res; + double cutinv2 = 1/cutmax/cutmax; + ilist = sim->ilist; + int nelements=pair->nelements; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double expr[jnum][kmax+12]; + int p = kmax; + int countmb=((mlength)*(mlength+1))>>1; + double *si = sim->s[i]; + // calculate interpolation expr, rinvs and dfc, for each neighbor + for (jj = 0; jj < jnum; jj++) { + if (Bij[jj]==false) {continue;} + jtype = tn[jj]; + if (atomtypes[1] != nelements && atomtypes[1] != jtype && atomtypes[2] != nelements && atomtypes[2] != jtype) { + expr[jj][0]=0; + continue; + } + delx=xn[jj]; + dely=yn[jj]; + delz=zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq>rc*rc) { + expr[jj][0]=0; + continue; + } + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + r1 = r1-trunc(r1); + double *p0 = &expcuttable[(m1-1)*kmax]; + double *p1 = &expcuttable[m1*kmax]; + double *p2 = &expcuttable[(m1+1)*kmax]; + double *p3 = &expcuttable[(m1+2)*kmax]; + for (kk=0;kks[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + double yprod = expr[jj][ai]; + double *y4 = &expr[jj][p]; + for (a2=0;a2s[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + double yprod = expr[jj][ai]; + double *y4 = &expr[jj][p]; + for (a2=0;a2s[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + double *y3 = &expr[jj][p+3]; + double *y4 = &expr[jj][p]; + ai = n; + yprod = expr[jj][ai]; + for (a2=0;a2s[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + double *y3 = &expr[jj][p+3]; + double *y4 = &expr[jj][p]; + ai = n; + yprod = expr[jj][ai]; + for (a2=0;a2s[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + double *y3 = &expr[jj][p+3]; + double *y4 = &expr[jj][p]; + ai = n; + yprod = expr[jj][ai]; + for (a2=0;a2sims[sid]; + int *type = sim->type; + int nelements = pair->nelements; + int res = pair->res; + double cutmax = pair->cutmax; + double cutinv2 = 1/cutmax/cutmax; + ilist = sim->ilist; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double expr[jnum][kmax]; + double y[jnum][3]; + double ri[jnum]; + double dfc[jnum]; + int kb = kmax; + int mb = mlength; + double Bijk[kb][mb]; + double c41[kmax]; + double c51[kmax]; + double c61[kmax]; + double ct[kmax]; + double *si = sim->s[i]; + for (jj = 0; jj < jnum; jj++) { + if (Bij[jj]==false) {continue;} + jtype = tn[jj]; + if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype) { + expr[jj][0]=0; + continue; + } + delx = xn[jj]; + dely = yn[jj]; + delz = zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq>rc*rc) { + expr[jj][0]=0; + continue; + } + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + if (!(m1>=1 && m1 <= res))pair->errorf(FLERR,"Neighbor list is invalid.");//usually results from nan somewhere. + r1 = r1-trunc(r1); + double *p0 = &expcuttable[(m1-1)*kmax]; + double *p1 = &expcuttable[m1*kmax]; + double *p2 = &expcuttable[(m1+1)*kmax]; + double *p3 = &expcuttable[(m1+2)*kmax]; + for (kk=0;kks[j]; + double spj = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + for (n = 0;ns[j]; + double spk = si[0]*sk[0]+si[1]*sk[1]+si[2]*sk[2]; + count = startingneuron; + double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]); + double c1 = ri[jj]*(y[kk][0]-dot*y[jj][0]); + double c2 = ri[jj]*(y[kk][1]-dot*y[jj][1]); + double c3 = ri[jj]*(y[kk][2]-dot*y[jj][2]); + double c10 = ri[kk]*(y[jj][0]-dot*y[kk][0]); + double c11 = ri[kk]*(y[jj][1]-dot*y[kk][1]); + double c12 = ri[kk]*(y[jj][2]-dot*y[kk][2]); + for (n=0;n); + void write_values(FILE *); + void init(int*,int); + void allocate(); + void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*); + void do3bodyfeatureset_doubleneighborloop(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int *); + void do3bodyfeatureset_singleneighborloop(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*); + void generate_exp_cut_table(); + void generate_coefficients(); + int get_length(); + + double *expcuttable; + double *dfctable; + double dr; + double *alpha_k; + double re; + int **coeff; + int **coeffx; + int **coeffy; + int **coeffz; + int kmax; + int mlength; + int **Mf; + + }; + + } +} + + +#endif +#endif /* FINGERPRINT_BOND_H_ */ diff --git a/src/USER-RANN/rann_fingerprint_bondspin.cpp b/src/USER-RANN/rann_fingerprint_bondspin.cpp new file mode 100644 index 0000000000..6aadbada86 --- /dev/null +++ b/src/USER-RANN/rann_fingerprint_bondspin.cpp @@ -0,0 +1,887 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ +#include "rann_fingerprint_bondspin.h" + + +using namespace LAMMPS_NS::RANN; + +Fingerprint_bondspin::Fingerprint_bondspin(PairRANN *_pair) : Fingerprint(_pair) +{ + n_body_type = 3; + dr = 0; + re = 0; + rc = 0; + alpha_k = new double [1]; + alpha_k[0] = -1; + kmax = 0; + mlength = 0; + id = -1; + style = "bondspin"; + atomtypes = new int [n_body_type]; + empty = true; + _pair->allscreen = false; + _pair->dospin = true; + spin = true; +} + +Fingerprint_bondspin::~Fingerprint_bondspin() { + delete [] alpha_k; + delete [] atomtypes; + delete [] expcuttable; + delete [] dfctable; + for (int i=0;i<(mlength*(mlength+1))>>1;i++) { + delete [] coeff[i]; + delete [] coeffx[i]; + delete [] coeffy[i]; + delete [] coeffz[i]; + delete [] Mf[i]; + } + delete [] coeff; + delete [] coeffx; + delete [] coeffy; + delete [] coeffz; + delete [] Mf; + delete [] rinvsqrttable; +} + +bool Fingerprint_bondspin::parse_values(std::string constant,std::vector line1) { + int nwords,l; + nwords=line1.size(); + if (constant.compare("re")==0) { + re = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("rc")==0) { + rc = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("alphak")==0) { + delete [] alpha_k; + alpha_k = new double [nwords]; + for (l=0;lerrorf(FLERR,"Undefined value for bond power"); + if (re!=0.0 && rc!=0.0 && alpha_k[0]!=-1 && dr!=0.0 && mlength!=0 && kmax!=0)return true; + return false; +} + +void Fingerprint_bondspin::write_values(FILE *fid) { + int i; + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:re:\n",style,id); + fprintf(fid,"%f\n",re); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:rc:\n",style,id); + fprintf(fid,"%f\n",rc); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:alphak:\n",style,id); + for (i=0;ielementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:dr:\n",style,id); + fprintf(fid,"%f\n",dr); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:k:\n",style,id); + fprintf(fid,"%d\n",kmax); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:m:\n",style,id); + fprintf(fid,"%d\n",mlength); +} + +void Fingerprint_bondspin::init(int *i,int _id) { + for (int j=0;jres; + double cutmax = pair->cutmax; + expcuttable = new double [(res+buf)*(kmax)]; + dfctable = new double [res+buf]; + for (m=0;m<(res+buf);m++) { + r1 = cutmax*cutmax*(double)(m)/(double)(res); + for (n=0;n<(kmax);n++) { + expcuttable[n+m*(kmax)] = exp(-alpha_k[n]/re*sqrt(r1))*cutofffunction(sqrt(r1),rc,dr); + } + if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) { + dfctable[m]=0; + } + else{ + dfctable[m]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4)); + } + } +} + +//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop. +void Fingerprint_bondspin::generate_coefficients() { //calculates multinomial coefficient for each term + int p,mb,mc; + int n,p1,i1; + mb = mlength; + mc=(mb*(mb+1))>>1; + coeff = new int *[mc]; + coeffx = new int *[mc]; + coeffy = new int *[mc]; + coeffz = new int *[mc]; + for (p=0;pfactorial(p)/pair->factorial(coeffx[p1][p])/pair->factorial(coeffy[p1][p])/pair->factorial(coeffz[p1][p]); + } + } +} + + +//Called by getproperties. Gets 3-body features and dfeatures +void Fingerprint_bondspin::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double * dspinx,double *dspiny,double *dspinz,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + int i; + int *ilist,*numneigh; + PairRANN::Simulation *sim = &pair->sims[sid]; + ilist = sim->ilist; + numneigh = sim->numneigh; + i = ilist[ii]; + //select the more efficient algorithm for this particular potential and environment. + if (jnum*2>(mlength+1)*mlength*20) { + do3bodyfeatureset_singleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,dspinx,dspiny,dspinz,ii,sid,xn,yn,zn,tn,jnum,jl); + } + else{ + do3bodyfeatureset_doubleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,dspinx,dspiny,dspinz,ii,sid,xn,yn,zn,tn,jnum,jl); + } +} + +//Called by do3bodyfeatureset. Algorithm for high neighbor numbers and small series of bond angle powers +void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double * dspinx,double *dspiny,double *dspinz,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) { + int i,j,jj,itype,jtype,kk,m,n,mcount,a,a1,a2,ai; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *ilist,*jlist,*numneigh,**firstneigh; + int count=0; + PairRANN::Simulation *sim = &pair->sims[sid]; + int *type = sim->type; + double cutmax = pair->cutmax; + int res = pair->res; + double cutinv2 = 1/cutmax/cutmax; + ilist = sim->ilist; + int nelements=pair->nelements; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double expr[jnum][kmax+12]; + int p = kmax; + int countmb=((mlength)*(mlength+1))>>1; + double *si = sim->s[i]; + // calculate interpolation expr, rinvs and dfc, for each neighbor + for (jj = 0; jj < jnum; jj++) { + jtype = tn[jj]; + if (atomtypes[1] != nelements && atomtypes[1] != jtype && atomtypes[2] != nelements && atomtypes[2] != jtype) { + expr[jj][0]=0; + continue; + } + delx = xn[jj]; + dely = yn[jj]; + delz = zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq>rc*rc) { + expr[jj][0]=0; + continue; + } + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + r1 = r1-trunc(r1); + double *p0 = &expcuttable[(m1-1)*kmax]; + double *p1 = &expcuttable[m1*kmax]; + double *p2 = &expcuttable[(m1+1)*kmax]; + double *p3 = &expcuttable[(m1+2)*kmax]; + for (kk=0;kks[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + double yprod = expr[jj][ai]; + double *y4 = &expr[jj][p]; + for (a2=0;a2s[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + double yprod = expr[jj][ai]; + double *y4 = &expr[jj][p]; + for (a2=0;a2s[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + double *y3 = &expr[jj][p+3]; + double *y4 = &expr[jj][p]; + ai = n; + yprod = expr[jj][ai]; + for (a2=0;a2s[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + double *y3 = &expr[jj][p+3]; + double *y4 = &expr[jj][p]; + ai = n; + yprod = expr[jj][ai]; + for (a2=0;a2s[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + double *y3 = &expr[jj][p+3]; + double *y4 = &expr[jj][p]; + ai = n; + yprod = expr[jj][ai]; + for (a2=0;a2sims[sid]; + double **x = sim->x; + int *type = sim->type; + int nelements = pair->nelements; + int res = pair->res; + double cutmax = pair->cutmax; + double cutinv2 = 1/cutmax/cutmax; + ilist = sim->ilist; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double expr[jnum][kmax]; + double y[jnum][3]; + double ri[jnum]; + double dfc[jnum]; + int kb = kmax; + int mb = mlength; + double c41[kmax]; + double c51[kmax]; + double c61[kmax]; + double ct[kmax]; + double *si = sim->s[i]; + for (jj = 0; jj < jnum; jj++) { + jtype = tn[jj]; + if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype) { + expr[jj][0]=0; + continue; + } + delx = xn[jj]; + dely = yn[jj]; + delz = zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq>rc*rc) { + expr[jj][0]=0; + continue; + } + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + if (!(m1>=1 && m1 <= res))pair->errorf(FLERR,"Neighbor list is invalid.");//usually results from nan somewhere. + r1 = r1-trunc(r1); + double *p0 = &expcuttable[(m1-1)*kmax]; + double *p1 = &expcuttable[m1*kmax]; + double *p2 = &expcuttable[(m1+1)*kmax]; + double *p3 = &expcuttable[(m1+2)*kmax]; + for (kk=0;kks[j]; + double spj = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + for (n = 0;ns[j]; + double spk = si[0]*sk[0]+si[1]*sk[1]+si[2]*sk[2]; + count = startingneuron; + double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]); + double c1 = 2*ri[jj]*(y[kk][0]-dot*y[jj][0]); + double c2 = 2*ri[jj]*(y[kk][1]-dot*y[jj][1]); + double c3 = 2*ri[jj]*(y[kk][2]-dot*y[jj][2]); + double c10 = 2*ri[kk]*(y[jj][0]-dot*y[kk][0]); + double c11 = 2*ri[kk]*(y[jj][1]-dot*y[kk][1]); + double c12 = 2*ri[kk]*(y[jj][2]-dot*y[kk][2]); + for (n=0;ns[j]; + double spk = si[0]*sk[0]+si[1]*sk[1]+si[2]*sk[2]; + count = startingneuron; + double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]); + double c1 = ri[jj]*(y[kk][0]-dot*y[jj][0]); + double c2 = ri[jj]*(y[kk][1]-dot*y[jj][1]); + double c3 = ri[jj]*(y[kk][2]-dot*y[jj][2]); + double c10 = ri[kk]*(y[jj][0]-dot*y[kk][0]); + double c11 = ri[kk]*(y[jj][1]-dot*y[kk][1]); + double c12 = ri[kk]*(y[jj][2]-dot*y[kk][2]); + for (n=0;n); + void write_values(FILE *); + void init(int*,int); + void allocate(); + virtual void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*);//spin + void do3bodyfeatureset_doubleneighborloop(double*,double*,double*,double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*); + void do3bodyfeatureset_singleneighborloop(double*,double*,double*,double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*); + void generate_exp_cut_table(); + void generate_coefficients(); + int get_length(); + + double *expcuttable; + double *dfctable; + double dr; + double *alpha_k; + double re; + int **coeff; + int **coeffx; + int **coeffy; + int **coeffz; + int kmax; + int mlength; + int **Mf; + + }; + + } +} + + +#endif +#endif /* FINGERPRINT_BOND_H_ */ diff --git a/src/USER-RANN/rann_fingerprint_radial.cpp b/src/USER-RANN/rann_fingerprint_radial.cpp new file mode 100644 index 0000000000..5df13aa094 --- /dev/null +++ b/src/USER-RANN/rann_fingerprint_radial.cpp @@ -0,0 +1,239 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ + +#include "rann_fingerprint_radial.h" + + + +using namespace LAMMPS_NS::RANN; + +Fingerprint_radial::Fingerprint_radial(PairRANN *_pair) : Fingerprint(_pair) +{ + n_body_type = 2; + dr = 0; + re = 0; + rc = 0; + alpha = new double [1]; + alpha[0] = -1; + nmax = 0; + omin = 0; + id = -1; + style = "radial"; + atomtypes = new int [n_body_type]; + empty = true; + fullydefined = false; + _pair->allscreen = false; +} + +Fingerprint_radial::~Fingerprint_radial() +{ + delete [] atomtypes; + delete [] radialtable; + delete [] alpha; + delete [] dfctable; + delete [] rinvsqrttable; +} + +bool Fingerprint_radial::parse_values(std::string constant,std::vector line1) { + int l; + int nwords=line1.size(); + if (constant.compare("re")==0) { + re = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("rc")==0) { + rc = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("alpha")==0) { + delete [] alpha; + alpha = new double [nwords]; + for (l=0;lerrorf(FLERR,"Undefined value for radial power"); + //code will run with default o=0 if o is never specified. All other values must be defined in potential file. + if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && nmax!=0)return true; + return false; +} + +void Fingerprint_radial::write_values(FILE *fid) { + int i; + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:re:\n",style,id); + fprintf(fid,"%f\n",re); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:rc:\n",style,id); + fprintf(fid,"%f\n",rc); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:alpha:\n",style,id); + for (i=0;i<(nmax-omin+1);i++) { + fprintf(fid,"%f ",alpha[i]); + } + fprintf(fid,"\n"); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:dr:\n",style,id); + fprintf(fid,"%f\n",dr); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:o:\n",style,id); + fprintf(fid,"%d\n",omin); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:n:\n",style,id); + fprintf(fid,"%d\n",nmax); +} + +//called after fingerprint is fully defined and tables can be computed. +void Fingerprint_radial::allocate() +{ + int k,m; + double r1; + int buf = 5; + int res = pair->res; + double cutmax = pair->cutmax; + radialtable = new double [(res+buf)*get_length()]; + dfctable = new double [res+buf]; + for (k=0;k<(res+buf);k++) { + r1 = cutmax*cutmax*(double)(k)/(double)(res); + for (m=0;m<=(nmax-omin);m++) { + radialtable[k*(nmax-omin+1)+m]=pow(sqrt(r1)/re,m+omin)*exp(-alpha[m]*sqrt(r1)/re)*cutofffunction(sqrt(r1),rc,dr); + } + if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) { + dfctable[k]=0; + } + else{ + dfctable[k]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4)); + } + } + generate_rinvssqrttable(); +} + +//called after fingerprint is declared for i-j type, but before its parameters are read. +void Fingerprint_radial::init(int *i,int _id) +{ + empty = false; + for (int j=0;jnelements; + int res = pair->res; + int i,j,jj,itype,jtype,l; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *ilist,*jlist,*numneigh,**firstneigh; + // + PairRANN::Simulation *sim = &pair->sims[sid]; + int count=0; +// double **x = sim->x; + int *type = sim->type; + ilist = sim->ilist; + double cutmax = pair->cutmax; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double cutinv2 = 1/cutmax/cutmax; + //loop over neighbors + for (jj = 0; jj < jnum; jj++) { + jtype =tn[jj]; + if (atomtypes[1] != nelements && atomtypes[1] != jtype)continue; + delx = xn[jj]; + dely = yn[jj]; + delz = zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq > rc*rc)continue; + count = startingneuron; + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + if (m1>res || m1<1) {pair->errorf(FLERR,"invalid neighbor radius!");} + if (radialtable[m1]==0) {continue;} + //cubic interpolation from tables + double *p1 = &radialtable[m1*(nmax-omin+1)]; + double *p2 = &radialtable[(m1+1)*(nmax-omin+1)]; + double *p3 = &radialtable[(m1+2)*(nmax-omin+1)]; + double *p0 = &radialtable[(m1-1)*(nmax-omin+1)]; + double *q = &dfctable[m1-1]; + double *rinvs = &rinvsqrttable[m1-1]; + r1 = r1-trunc(r1); + double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0]))); + double ri = rinvs[1] + 0.5 * r1*(rinvs[2] - rinvs[0] + r1*(2.0*rinvs[0] - 5.0*rinvs[1] + 4.0*rinvs[2] - rinvs[3] + r1*(3.0*(rinvs[1] - rinvs[2]) + rinvs[3] - rinvs[0]))); + for (l=0;l<=(nmax-omin);l++) { + double rt = p1[l]+0.5*r1*(p2[l]-p0[l]+r1*(2.0*p0[l]-5.0*p1[l]+4.0*p2[l]-p3[l]+r1*(3.0*(p1[l]-p2[l])+p3[l]-p0[l]))); + features[count]+=rt; + rt *= (l+omin)/rsq+(-alpha[l]/re+dfc)*ri; + //update neighbor's features + dfeaturesx[jj*f+count]+=rt*delx; + dfeaturesy[jj*f+count]+=rt*dely; + dfeaturesz[jj*f+count]+=rt*delz; + //update atom's features + dfeaturesx[jnum*f+count]-=rt*delx; + dfeaturesy[jnum*f+count]-=rt*dely; + dfeaturesz[jnum*f+count]-=rt*delz; + count++; + } + } +} + +int Fingerprint_radial::get_length() +{ + return nmax-omin+1; +} diff --git a/src/USER-RANN/rann_fingerprint_radial.h b/src/USER-RANN/rann_fingerprint_radial.h new file mode 100644 index 0000000000..5f9876f493 --- /dev/null +++ b/src/USER-RANN/rann_fingerprint_radial.h @@ -0,0 +1,70 @@ + +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ +#ifdef FINGERPRINT_CLASS + +FingerprintStyle(radial,Fingerprint_radial) + +#else + +#ifndef LMP_RANN_FINGERPRINT_RADIAL_H +#define LMP_RANN_FINGERPRINT_RADIAL_H + +#include "rann_fingerprint.h" + +namespace LAMMPS_NS { + namespace RANN { + + class Fingerprint_radial : public Fingerprint { + public: + Fingerprint_radial(PairRANN *); + ~Fingerprint_radial(); + bool parse_values(std::string,std::vector); + void write_values(FILE *); + void init(int*,int); + void allocate(); + void compute_fingerprint(double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*); + int get_length(); + + double *radialtable; + double *dfctable; + double dr; + double *alpha; + double re; + int nmax;//highest term + int omin;//lowest term + + }; + + } +} + +#endif +#endif /* FINGERPRINT_RADIAL_H_ */ diff --git a/src/USER-RANN/rann_fingerprint_radialscreened.cpp b/src/USER-RANN/rann_fingerprint_radialscreened.cpp new file mode 100644 index 0000000000..79333a4e3a --- /dev/null +++ b/src/USER-RANN/rann_fingerprint_radialscreened.cpp @@ -0,0 +1,253 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ + +#include "rann_fingerprint_radialscreened.h" + + + +using namespace LAMMPS_NS::RANN; + +Fingerprint_radialscreened::Fingerprint_radialscreened(PairRANN *_pair) : Fingerprint(_pair) +{ + n_body_type = 2; + dr = 0; + re = 0; + rc = 0; + alpha = new double [1]; + alpha[0] = -1; + nmax = 0; + omin = 0; + id = -1; + style = "radialscreened"; + atomtypes = new int [n_body_type]; + empty = true; + fullydefined = false; + _pair->doscreen = true; + screen = true; +} + +Fingerprint_radialscreened::~Fingerprint_radialscreened() +{ + delete [] atomtypes; + delete [] radialtable; + delete [] alpha; + delete [] dfctable; + delete [] rinvsqrttable; +} + +bool Fingerprint_radialscreened::parse_values(std::string constant,std::vector line1) { + int l; + int nwords=line1.size(); + if (constant.compare("re")==0) { + re = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("rc")==0) { + rc = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("alpha")==0) { + delete [] alpha; + alpha = new double [nwords]; + for (l=0;lerrorf(FLERR,"Undefined value for radial power"); + //code will run with default o=0 if o is never specified. All other values must be defined in potential file. + if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && nmax!=0)return true; + return false; +} + +void Fingerprint_radialscreened::write_values(FILE *fid) { + int i; + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:re:\n",style,id); + fprintf(fid,"%f\n",re); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:rc:\n",style,id); + fprintf(fid,"%f\n",rc); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:alpha:\n",style,id); + for (i=0;i<(nmax-omin+1);i++) { + fprintf(fid,"%f ",alpha[i]); + } + fprintf(fid,"\n"); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:dr:\n",style,id); + fprintf(fid,"%f\n",dr); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:o:\n",style,id); + fprintf(fid,"%d\n",omin); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:n:\n",style,id); + fprintf(fid,"%d\n",nmax); +} + +//called after fingerprint is fully defined and tables can be computed. +void Fingerprint_radialscreened::allocate() +{ + int k,m; + double r1; + int buf = 5; + int res = pair->res; + double cutmax = pair->cutmax; + radialtable = new double [(res+buf)*get_length()]; + dfctable = new double [res+buf]; + for (k=0;k<(res+buf);k++) { + r1 = cutmax*cutmax*(double)(k)/(double)(res); + for (m=0;m<=(nmax-omin);m++) { + radialtable[k*(nmax-omin+1)+m]=pow(sqrt(r1)/re,m+omin)*exp(-alpha[m]*sqrt(r1)/re)*cutofffunction(sqrt(r1),rc,dr); + } + if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) { + dfctable[k]=0; + } + else{ + dfctable[k]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4)); + } + } + generate_rinvssqrttable(); +} + +//called after fingerprint is declared for i-j type, but before its parameters are read. +void Fingerprint_radialscreened::init(int *i,int _id) +{ + empty = false; + for (int j=0;jnelements; + int res = pair->res; + int i,j,jj,itype,jtype,l,kk; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *ilist,*jlist,*numneigh,**firstneigh; + // + PairRANN::Simulation *sim = &pair->sims[sid]; + int count=0; + double **x = sim->x; + int *type = sim->type; + ilist = sim->ilist; + double cutmax = pair->cutmax; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double cutinv2 = 1/cutmax/cutmax; + //loop over neighbors + for (jj = 0; jj < jnum; jj++) { + if (Bij[jj]==false) {continue;} + jtype = tn[jj]; + if (atomtypes[1] != nelements && atomtypes[1] != jtype)continue; + delx = xn[jj]; + dely = yn[jj]; + delz = zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq > rc*rc)continue; + count = startingneuron; + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + if (m1>res || m1<1) {pair->errorf(FLERR,"invalid neighbor radius!");} + if (radialtable[m1]==0) {continue;} + //cubic interpolation from tables + double *p1 = &radialtable[m1*(nmax-omin+1)]; + double *p2 = &radialtable[(m1+1)*(nmax-omin+1)]; + double *p3 = &radialtable[(m1+2)*(nmax-omin+1)]; + double *p0 = &radialtable[(m1-1)*(nmax-omin+1)]; + double *q = &dfctable[m1-1]; + double *rinvs = &rinvsqrttable[m1-1]; + r1 = r1-trunc(r1); + double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0]))); + double ri = rinvs[1] + 0.5 * r1*(rinvs[2] - rinvs[0] + r1*(2.0*rinvs[0] - 5.0*rinvs[1] + 4.0*rinvs[2] - rinvs[3] + r1*(3.0*(rinvs[1] - rinvs[2]) + rinvs[3] - rinvs[0]))); + for (l=0;l<=(nmax-omin);l++) { + double rt = Sik[jj]*(p1[l]+0.5*r1*(p2[l]-p0[l]+r1*(2.0*p0[l]-5.0*p1[l]+4.0*p2[l]-p3[l]+r1*(3.0*(p1[l]-p2[l])+p3[l]-p0[l])))); + features[count]+=rt; + double rt1 = rt*((l+omin)/rsq+(-alpha[l]/re+dfc)*ri); + //update neighbor's features + dfeaturesx[jj*f+count]+=rt1*delx+rt*dSikx[jj]; + dfeaturesy[jj*f+count]+=rt1*dely+rt*dSiky[jj]; + dfeaturesz[jj*f+count]+=rt1*delz+rt*dSikz[jj]; + for (kk=0;kk); + void write_values(FILE *); + void init(int*,int); + void allocate(); + void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*); + int get_length(); + + double *radialtable; + double *dfctable; + double dr; + double *alpha; + double re; + int nmax;//highest term + int omin;//lowest term + + }; +} + +} + +#endif + + + +#endif /* FINGERPRINT_RADIALSCREENED_H_ */ diff --git a/src/USER-RANN/rann_fingerprint_radialscreenedspin.cpp b/src/USER-RANN/rann_fingerprint_radialscreenedspin.cpp new file mode 100644 index 0000000000..f594e92077 --- /dev/null +++ b/src/USER-RANN/rann_fingerprint_radialscreenedspin.cpp @@ -0,0 +1,264 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ +#include "rann_fingerprint_radialscreenedspin.h" + + + +using namespace LAMMPS_NS::RANN; + +Fingerprint_radialscreenedspin::Fingerprint_radialscreenedspin(PairRANN *_pair) : Fingerprint(_pair) +{ + n_body_type = 2; + dr = 0; + re = 0; + rc = 0; + alpha = new double [1]; + alpha[0] = -1; + nmax = 0; + omin = 0; + id = -1; + style = "radialscreenedspin"; + atomtypes = new int [n_body_type]; + empty = true; + fullydefined = false; + _pair->doscreen = true; + screen = true; + _pair->dospin = true; + spin = true; +} + +Fingerprint_radialscreenedspin::~Fingerprint_radialscreenedspin() +{ + delete [] atomtypes; + delete [] radialtable; + delete [] alpha; + delete [] dfctable; + delete [] rinvsqrttable; +} + +bool Fingerprint_radialscreenedspin::parse_values(std::string constant,std::vector line1) { + int l; + int nwords=line1.size(); + if (constant.compare("re")==0) { + re = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("rc")==0) { + rc = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("alpha")==0) { + delete [] alpha; + alpha = new double [nwords]; + for (l=0;lerrorf(FLERR,"Undefined value for radial power"); + //code will run with default o=0 if o is never specified. All other values must be defined in potential file. + if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && nmax!=0)return true; + return false; +} + +void Fingerprint_radialscreenedspin::write_values(FILE *fid) { + int i; + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:re:\n",style,id); + fprintf(fid,"%f\n",re); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:rc:\n",style,id); + fprintf(fid,"%f\n",rc); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:alpha:\n",style,id); + for (i=0;i<(nmax-omin+1);i++) { + fprintf(fid,"%f ",alpha[i]); + } + fprintf(fid,"\n"); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:dr:\n",style,id); + fprintf(fid,"%f\n",dr); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:o:\n",style,id); + fprintf(fid,"%d\n",omin); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:n:\n",style,id); + fprintf(fid,"%d\n",nmax); +} + +//called after fingerprint is fully defined and tables can be computed. +void Fingerprint_radialscreenedspin::allocate() +{ + int k,m; + double r1; + int buf = 5; + int res = pair->res; + double cutmax = pair->cutmax; + radialtable = new double [(res+buf)*get_length()]; + dfctable = new double [res+buf]; + for (k=0;k<(res+buf);k++) { + r1 = cutmax*cutmax*(double)(k)/(double)(res); + for (m=0;m<=(nmax-omin);m++) { + radialtable[k*(nmax-omin+1)+m]=pow(sqrt(r1)/re,m+omin)*exp(-alpha[m]*sqrt(r1)/re)*cutofffunction(sqrt(r1),rc,dr); + } + if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) { + dfctable[k]=0; + } + else{ + dfctable[k]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4)); + } + } + generate_rinvssqrttable(); +} + +//called after fingerprint is declared for i-j type, but before its parameters are read. +void Fingerprint_radialscreenedspin::init(int *i,int _id) +{ + empty = false; + for (int j=0;jnelements; + int res = pair->res; + int i,j,jj,itype,jtype,l,kk; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *ilist,*jlist,*numneigh,**firstneigh; + PairRANN::Simulation *sim = &pair->sims[sid]; + int count=0; + double **x = sim->x; + int *type = sim->type; + ilist = sim->ilist; + double cutmax = pair->cutmax; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double cutinv2 = 1/cutmax/cutmax; + double *si = sim->s[i]; + //loop over neighbors + for (jj = 0; jj < jnum; jj++) { + if (Bij[jj]==false) {continue;} + jtype = tn[jj]; + if (atomtypes[1] != nelements && atomtypes[1] != jtype)continue; + delx = xn[jj]; + dely = yn[jj]; + delz = zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq > rc*rc)continue; + count = startingneuron; + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + if (m1>res || m1<1) {pair->errorf(FLERR,"invalid neighbor radius!");} + if (radialtable[m1]==0) {continue;} + j=jl[jj]; + double *sj = sim->s[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + //cubic interpolation from tables + double *p1 = &radialtable[m1*(nmax-omin+1)]; + double *p2 = &radialtable[(m1+1)*(nmax-omin+1)]; + double *p3 = &radialtable[(m1+2)*(nmax-omin+1)]; + double *p0 = &radialtable[(m1-1)*(nmax-omin+1)]; + double *q = &dfctable[m1-1]; + double *rinvs = &rinvsqrttable[m1-1]; + r1 = r1-trunc(r1); + double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0]))); + double ri = rinvs[1] + 0.5 * r1*(rinvs[2] - rinvs[0] + r1*(2.0*rinvs[0] - 5.0*rinvs[1] + 4.0*rinvs[2] - rinvs[3] + r1*(3.0*(rinvs[1] - rinvs[2]) + rinvs[3] - rinvs[0]))); + for (l=0;l<=(nmax-omin);l++) { + double rt = Sik[jj]*(p1[l]+0.5*r1*(p2[l]-p0[l]+r1*(2.0*p0[l]-5.0*p1[l]+4.0*p2[l]-p3[l]+r1*(3.0*(p1[l]-p2[l])+p3[l]-p0[l])))); + //update neighbor's features + dspinx[jj*f+count]+=rt*si[0]; + dspiny[jj*f+count]+=rt*si[1]; + dspinz[jj*f+count]+=rt*si[2]; + dspinx[jnum*f+count]+=rt*sj[0]; + dspiny[jnum*f+count]+=rt*sj[1]; + dspinz[jnum*f+count]+=rt*sj[2]; + rt *= sp; + features[count]+=rt; + double rt1 = rt*((l+omin)/rsq+(-alpha[l]/re+dfc)*ri); + dfeaturesx[jj*f+count]+=rt1*delx+rt*dSikx[jj]; + dfeaturesy[jj*f+count]+=rt1*dely+rt*dSiky[jj]; + dfeaturesz[jj*f+count]+=rt1*delz+rt*dSikz[jj]; + for (kk=0;kk); + void write_values(FILE *); + void init(int*,int); + void allocate(); + virtual void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);//spin,screen + int get_length(); + + double *radialtable; + double *dfctable; + double dr; + double *alpha; + double re; + int nmax;//highest term + int omin;//lowest term + + }; + } + +} + +#endif + + + +#endif /* FINGERPRINT_RADIALSCREENED_H_ */ diff --git a/src/USER-RANN/rann_fingerprint_radialspin.cpp b/src/USER-RANN/rann_fingerprint_radialspin.cpp new file mode 100644 index 0000000000..9d61fc08ed --- /dev/null +++ b/src/USER-RANN/rann_fingerprint_radialspin.cpp @@ -0,0 +1,252 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ +#include "rann_fingerprint_radialspin.h" + + + +using namespace LAMMPS_NS::RANN; + +Fingerprint_radialspin::Fingerprint_radialspin(PairRANN *_pair) : Fingerprint(_pair) +{ + n_body_type = 2; + dr = 0; + re = 0; + rc = 0; + alpha = new double [1]; + alpha[0] = -1; + nmax = 0; + omin = 0; + id = -1; + style = "radialspin"; + atomtypes = new int [n_body_type]; + empty = true; + fullydefined = false; + _pair->allscreen = false; + _pair->dospin = true; + spin = true; +} + +Fingerprint_radialspin::~Fingerprint_radialspin() +{ + delete [] atomtypes; + delete [] radialtable; + delete [] alpha; + delete [] dfctable; + delete [] rinvsqrttable; +} + +bool Fingerprint_radialspin::parse_values(std::string constant,std::vector line1) { + int l; + int nwords=line1.size(); + if (constant.compare("re")==0) { + re = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("rc")==0) { + rc = strtod(line1[0].c_str(),NULL); + } + else if (constant.compare("alpha")==0) { + delete [] alpha; + alpha = new double [nwords]; + for (l=0;lerrorf(FLERR,"Undefined value for radial power"); + //code will run with default o=0 if o is never specified. All other values must be defined in potential file. + if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && nmax!=0)return true; + return false; +} + +void Fingerprint_radialspin::write_values(FILE *fid) { + int i; + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:re:\n",style,id); + fprintf(fid,"%f\n",re); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:rc:\n",style,id); + fprintf(fid,"%f\n",rc); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:alpha:\n",style,id); + for (i=0;i<(nmax-omin+1);i++) { + fprintf(fid,"%f ",alpha[i]); + } + fprintf(fid,"\n"); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:dr:\n",style,id); + fprintf(fid,"%f\n",dr); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:o:\n",style,id); + fprintf(fid,"%d\n",omin); + fprintf(fid,"fingerprintconstants:"); + fprintf(fid,"%s",pair->elementsp[atomtypes[0]]); + for (i=1;ielementsp[atomtypes[i]]); + } + fprintf(fid,":%s_%d:n:\n",style,id); + fprintf(fid,"%d\n",nmax); +} + +//called after fingerprint is fully defined and tables can be computed. +void Fingerprint_radialspin::allocate() +{ + int k,m; + double r1; + int buf = 5; + int res = pair->res; + double cutmax = pair->cutmax; + radialtable = new double [(res+buf)*get_length()]; + dfctable = new double [res+buf]; + for (k=0;k<(res+buf);k++) { + r1 = cutmax*cutmax*(double)(k)/(double)(res); + for (m=0;m<=(nmax-omin);m++) { + radialtable[k*(nmax-omin+1)+m]=pow(sqrt(r1)/re,m+omin)*exp(-alpha[m]*sqrt(r1)/re)*cutofffunction(sqrt(r1),rc,dr); + } + if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) { + dfctable[k]=0; + } + else{ + dfctable[k]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4)); + } + } + generate_rinvssqrttable(); +} + +//called after fingerprint is declared for i-j type, but before its parameters are read. +void Fingerprint_radialspin::init(int *i,int _id) +{ + empty = false; + for (int j=0;jnelements; + int res = pair->res; + int i,j,jj,itype,jtype,l; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *ilist,*jlist,*numneigh,**firstneigh; + // + PairRANN::Simulation *sim = &pair->sims[sid]; + int count=0; + int *type = sim->type; + ilist = sim->ilist; + double cutmax = pair->cutmax; + i = ilist[ii]; + itype = pair->map[type[i]]; + int f = pair->net[itype].dimensions[0]; + double cutinv2 = 1/cutmax/cutmax; + double *si = sim->s[i]; + firstneigh = sim->firstneigh; + jlist = firstneigh[i]; + //loop over neighbors + for (jj = 0; jj < jnum; jj++) { + j = jl[jj]; + jtype =tn[jj]; + if (atomtypes[1] != nelements && atomtypes[1] != jtype)continue; + delx = xn[jj]; + dely = yn[jj]; + delz = zn[jj]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq > rc*rc)continue; + count = startingneuron; + double r1 = (rsq*((double)res)*cutinv2); + int m1 = (int)r1; + if (m1>res || m1<1) {pair->errorf(FLERR,"invalid neighbor radius!");} + if (radialtable[m1]==0) {continue;} + double *sj = sim->s[j]; + double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2]; + //cubic interpolation from tables + double *p1 = &radialtable[m1*(nmax-omin+1)]; + double *p2 = &radialtable[(m1+1)*(nmax-omin+1)]; + double *p3 = &radialtable[(m1+2)*(nmax-omin+1)]; + double *p0 = &radialtable[(m1-1)*(nmax-omin+1)]; + double *q = &dfctable[m1-1]; + double *rinvs = &rinvsqrttable[m1-1]; + r1 = r1-trunc(r1); + double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0]))); + double ri = rinvs[1] + 0.5 * r1*(rinvs[2] - rinvs[0] + r1*(2.0*rinvs[0] - 5.0*rinvs[1] + 4.0*rinvs[2] - rinvs[3] + r1*(3.0*(rinvs[1] - rinvs[2]) + rinvs[3] - rinvs[0]))); + for (l=0;l<=(nmax-omin);l++) { + double rt = p1[l]+0.5*r1*(p2[l]-p0[l]+r1*(2.0*p0[l]-5.0*p1[l]+4.0*p2[l]-p3[l]+r1*(3.0*(p1[l]-p2[l])+p3[l]-p0[l]))); + dspinx[jj*f+count]+=rt*si[0]; + dspiny[jj*f+count]+=rt*si[1]; + dspinz[jj*f+count]+=rt*si[2]; + dspinx[jnum*f+count]+=rt*sj[0]; + dspiny[jnum*f+count]+=rt*sj[1]; + dspinz[jnum*f+count]+=rt*sj[2]; + rt *= sp; + features[count]+=rt; + rt *= (l+omin)/rsq+(-alpha[l]/re+dfc)*ri; + //update neighbor's features + dfeaturesx[jj*f+count]+=rt*delx; + dfeaturesy[jj*f+count]+=rt*dely; + dfeaturesz[jj*f+count]+=rt*delz; + //update atom's features + dfeaturesx[jnum*f+count]-=rt*delx; + dfeaturesy[jnum*f+count]-=rt*dely; + dfeaturesz[jnum*f+count]-=rt*delz; + count++; + } + } +} + +int Fingerprint_radialspin::get_length() +{ + return nmax-omin+1; +} diff --git a/src/USER-RANN/rann_fingerprint_radialspin.h b/src/USER-RANN/rann_fingerprint_radialspin.h new file mode 100644 index 0000000000..6abddd910d --- /dev/null +++ b/src/USER-RANN/rann_fingerprint_radialspin.h @@ -0,0 +1,69 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu + Doyl Dickel (MSU) doyl@me.msstate.edu + ----------------------------------------------------------------------*/ +/* +“The research described and the resulting data presented herein, unless +otherwise noted, was funded under PE 0602784A, Project T53 "Military +Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095, +managed by the U.S. Army Combat Capabilities Development Command (CCDC) and +the Engineer Research and Development Center (ERDC). The work described in +this document was conducted at CAVS, MSU. Permission was granted by ERDC +to publish this information. Any opinions, findings and conclusions or +recommendations expressed in this material are those of the author(s) and +do not necessarily reflect the views of the United States Army.​” + +DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918 + */ + +#ifdef FINGERPRINT_CLASS + +FingerprintStyle(radialspin,Fingerprint_radialspin) + +#else + +#ifndef LMP_RANN_FINGERPRINT_RADIALSPIN_H +#define LMP_RANN_FINGERPRINT_RADIALSPIN_H + +#include "rann_fingerprint.h" + +namespace LAMMPS_NS { + namespace RANN { + class Fingerprint_radialspin : public Fingerprint { + public: + Fingerprint_radialspin(PairRANN *); + ~Fingerprint_radialspin(); + bool parse_values(std::string,std::vector); + void write_values(FILE *); + void init(int*,int); + void allocate(); + void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*); + int get_length(); + + double *radialtable; + double *dfctable; + double dr; + double *alpha; + double re; + int nmax;//highest term + int omin;//lowest term + + }; + } + +} + +#endif +#endif /* FINGERPRINT_RADIAL_H_ */ diff --git a/src/USER-RANN/rann_list_activation.h b/src/USER-RANN/rann_list_activation.h new file mode 100644 index 0000000000..a140b92c65 --- /dev/null +++ b/src/USER-RANN/rann_list_activation.h @@ -0,0 +1,2 @@ +#include "rann_activation_linear.h" +#include "rann_activation_sig_i.h" diff --git a/src/USER-RANN/rann_list_fingerprint.h b/src/USER-RANN/rann_list_fingerprint.h new file mode 100644 index 0000000000..005c7f1b04 --- /dev/null +++ b/src/USER-RANN/rann_list_fingerprint.h @@ -0,0 +1,8 @@ +#include "rann_fingerprint_bond.h" +#include "rann_fingerprint_bondscreened.h" +#include "rann_fingerprint_bondscreenedspin.h" +#include "rann_fingerprint_bondspin.h" +#include "rann_fingerprint_radial.h" +#include "rann_fingerprint_radialscreened.h" +#include "rann_fingerprint_radialscreenedspin.h" +#include "rann_fingerprint_radialspin.h"