Initial inclusion of rann potential
This commit is contained in:
kipbarrett
44
src/activation.cpp
Normal file
44
src/activation.cpp
Normal file
@ -0,0 +1,44 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
|
||||
#include "activation.h"
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
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;
|
||||
}
|
||||
39
src/activation.h
Normal file
39
src/activation.h
Normal file
@ -0,0 +1,39 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#ifndef ACTIVATION_H_
|
||||
#define ACTIVATION_H_
|
||||
|
||||
#include "pair_rann.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
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;
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
|
||||
#endif /* ACTIVATION_H_ */
|
||||
42
src/activation_linear.cpp
Normal file
42
src/activation_linear.cpp
Normal file
@ -0,0 +1,42 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#include <math.h>
|
||||
#include "activation_linear.h"
|
||||
#include "activation.h"
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
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;
|
||||
}
|
||||
42
src/activation_linear.h
Normal file
42
src/activation_linear.h
Normal file
@ -0,0 +1,42 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
*/
|
||||
|
||||
#ifdef ACTIVATION_CLASS
|
||||
|
||||
ActivationStyle(linear,Activation_linear)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef ACTIVATION_LINEAR_H_
|
||||
#define ACTIVATION_LINEAR_H_
|
||||
|
||||
#include "activation.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Activation_linear : public Activation {
|
||||
public:
|
||||
Activation_linear(class PairRANN *);
|
||||
double activation_function(double);
|
||||
double dactivation_function(double);
|
||||
double ddactivation_function(double);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif /* ACTIVATION_LINEAR_H_ */
|
||||
41
src/activation_sigI.cpp
Normal file
41
src/activation_sigI.cpp
Normal file
@ -0,0 +1,41 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#include <math.h>
|
||||
#include "activation_sigI.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
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);;
|
||||
}
|
||||
42
src/activation_sigI.h
Normal file
42
src/activation_sigI.h
Normal file
@ -0,0 +1,42 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
*/
|
||||
|
||||
#ifdef ACTIVATION_CLASS
|
||||
|
||||
ActivationStyle(sigI,Activation_sigI)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef ACTIVATION_SIGI_H_
|
||||
#define ACTIVATION_SIGI_H_
|
||||
|
||||
#include "activation.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Activation_sigI : public Activation {
|
||||
public:
|
||||
Activation_sigI(class PairRANN *);
|
||||
double activation_function(double);
|
||||
double dactivation_function(double);
|
||||
double ddactivation_function(double);
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
#endif /* ACTIVATION_SIGI_H_ */
|
||||
106
src/fingerprint.cpp
Normal file
106
src/fingerprint.cpp
Normal file
@ -0,0 +1,106 @@
|
||||
/* -*- 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
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#include "fingerprint.h"
|
||||
#include <math.h>
|
||||
#include <cmath>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <iostream>
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
Fingerprint::Fingerprint(PairRANN *pair)
|
||||
{
|
||||
spin = false;
|
||||
screen = false;
|
||||
empty = true;
|
||||
fullydefined = false;
|
||||
n_body_type = 0;
|
||||
style = "empty";
|
||||
this->pair = pair;
|
||||
}
|
||||
|
||||
Fingerprint::~Fingerprint(){
|
||||
|
||||
}
|
||||
|
||||
bool Fingerprint::parse_values(char *word,char *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 = pow(1-pow(1-(rc-r)/dr,4.0),2.0);
|
||||
}
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
58
src/fingerprint.h
Normal file
58
src/fingerprint.h
Normal file
@ -0,0 +1,58 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#ifndef FINGERPRINT_H_
|
||||
#define FINGERPRINT_H_
|
||||
|
||||
#include "pair_rann.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Fingerprint {
|
||||
public:
|
||||
Fingerprint(PairRANN *);
|
||||
virtual ~Fingerprint();
|
||||
virtual bool parse_values(char*,char*);
|
||||
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 /* FINGERPRINT_H_ */
|
||||
912
src/fingerprint_bond.cpp
Normal file
912
src/fingerprint_bond.cpp
Normal file
@ -0,0 +1,912 @@
|
||||
/* -*- 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@cavs.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
|
||||
#include "fingerprint_bond.h"
|
||||
#include "fingerprint.h"
|
||||
#include <math.h>
|
||||
#include <cmath>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <iostream>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
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;
|
||||
k = 0;
|
||||
m = 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<(m*(m+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(char * constant, char * line1){
|
||||
char **words=new char *[MAXLINE];
|
||||
int nwords,l;
|
||||
nwords=0;
|
||||
words[nwords++] = strtok(line1,": ,\t\n");
|
||||
while ((words[nwords++] = strtok(NULL,": ,\t\n"))) continue;
|
||||
nwords -= 1;
|
||||
if (strcmp(constant,"re")==0){
|
||||
re = strtod(words[0],NULL);
|
||||
}
|
||||
else if (strcmp(constant,"rc")==0){
|
||||
rc = strtod(words[0],NULL);
|
||||
}
|
||||
else if (strcmp(constant,"alphak")==0){
|
||||
delete [] alpha_k;
|
||||
alpha_k = new double [nwords];
|
||||
for (l=0;l<nwords;l++){
|
||||
alpha_k[l]=strtod(words[l],NULL);
|
||||
}
|
||||
}
|
||||
else if (strcmp(constant,"dr")==0){
|
||||
dr = strtod(words[0],NULL);
|
||||
}
|
||||
else if (strcmp(constant,"k")==0){
|
||||
k = strtol(words[0],NULL,10);
|
||||
}
|
||||
else if (strcmp(constant,"m")==0){
|
||||
m = strtol(words[0],NULL,10);
|
||||
}
|
||||
else pair->errorf("Undefined value for bond power");
|
||||
delete [] words;
|
||||
if (re!=0.0 && rc!=0.0 && alpha_k[0]!=-1 && dr!=0.0 && m!=0 && k!=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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:alphak:\n",style,id);
|
||||
for (i=0;i<k;i++){
|
||||
fprintf(fid,"%f ",alpha_k[i]);
|
||||
}
|
||||
fprintf(fid,"\n");
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:k:\n",style,id);
|
||||
fprintf(fid,"%d\n",k);
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:m:\n",style,id);
|
||||
fprintf(fid,"%d\n",m);
|
||||
}
|
||||
|
||||
void Fingerprint_bond::init(int *i,int id){
|
||||
for (int j=0;j<n_body_type;j++){atomtypes[j] = i[j];}
|
||||
re = 0;
|
||||
rc = 0;
|
||||
m = 0;
|
||||
this->k = 0;
|
||||
alpha_k = new double [1];
|
||||
alpha_k[0]=-1;
|
||||
empty = false;
|
||||
this->id = id;
|
||||
}
|
||||
|
||||
//number of neurons defined by this fingerprint
|
||||
int Fingerprint_bond::get_length(){
|
||||
return m*k;
|
||||
}
|
||||
|
||||
void Fingerprint_bond::allocate(){
|
||||
generate_exp_cut_table();
|
||||
generate_coefficients();
|
||||
generate_rinvssqrttable();
|
||||
}
|
||||
|
||||
//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop and do3bodyfeatureset_doubleneighborloop.
|
||||
void Fingerprint_bond::generate_exp_cut_table(){
|
||||
int m,n;
|
||||
double r1;
|
||||
int buf = 5;
|
||||
int res = pair->res;
|
||||
double cutmax = pair->cutmax;
|
||||
expcuttable = new double [(res+buf)*(this->k)];
|
||||
dfctable = new double [res+buf];
|
||||
for (m=0;m<(res+buf);m++){
|
||||
r1 = cutmax*cutmax*(double)(m)/(double)(res);
|
||||
for (n=0;n<(this->k);n++){
|
||||
expcuttable[n+m*(this->k)] = 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 = this->m;
|
||||
mc=(mb*(mb+1))>>1;
|
||||
coeff = new int *[mc];
|
||||
coeffx = new int *[mc];
|
||||
coeffy = new int *[mc];
|
||||
coeffz = new int *[mc];
|
||||
for (p=0;p<mc;p++){
|
||||
coeff[p]=new int [mb];
|
||||
coeffx[p]=new int [mb];
|
||||
coeffy[p]=new int [mb];
|
||||
coeffz[p]=new int [mb];
|
||||
}
|
||||
Mf = new int*[mc];
|
||||
int *M = new int[this->m+1];
|
||||
for (p=0;p<this->m+1;p++){
|
||||
M[p]=0;
|
||||
}
|
||||
for (p1=0;p1<mc;p1++){
|
||||
Mf[p1] = new int[m+1];
|
||||
for (p=0;p<m+1;p++){
|
||||
Mf[p1][p]=0;
|
||||
}
|
||||
}
|
||||
M[0] = 2;
|
||||
Mf[0][0] = 2;
|
||||
n = 1;
|
||||
int m1 = 1;
|
||||
bool go = true;
|
||||
bool broke = false;
|
||||
while (go){
|
||||
broke = false;
|
||||
for (i1=0;i1<m-1;i1++){
|
||||
if (M[i1+1] == 0){
|
||||
M[i1+1]=M[i1+1]+1;
|
||||
for (p1=0;p1<m+1;p1++){
|
||||
Mf[n][p1] = M[p1];
|
||||
}
|
||||
n = n+1;
|
||||
broke = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (m1<m && !broke){
|
||||
M[m1]=M[m1]+1;
|
||||
for (p1=m1+1;p1<m+1;p1++){
|
||||
M[p1]=0;
|
||||
}
|
||||
for (p1=0;p1<m+1;p1++){
|
||||
Mf[n][p1]=M[p1];
|
||||
}
|
||||
n=n+1;
|
||||
broke = true;
|
||||
m1 = m1+1;
|
||||
}
|
||||
if (!broke){
|
||||
go = false;
|
||||
}
|
||||
}
|
||||
for (p=0;p<mb;p++){
|
||||
for (p1=0;p1<mc;p1++){
|
||||
if (p==0){
|
||||
coeffx[p1][p]=0;
|
||||
coeffy[p1][p]=0;
|
||||
coeffz[p1][p]=0;
|
||||
}
|
||||
else{
|
||||
coeffx[p1][p]=coeffx[p1][p-1];
|
||||
if (Mf[p1][p]==0){
|
||||
coeffx[p1][p]++;
|
||||
}
|
||||
coeffy[p1][p]=coeffy[p1][p-1];
|
||||
if (Mf[p1][p]==1){
|
||||
coeffy[p1][p]++;
|
||||
}
|
||||
coeffz[p1][p]=coeffz[p1][p-1];
|
||||
if (Mf[p1][p]==2){
|
||||
coeffz[p1][p]++;
|
||||
}
|
||||
}
|
||||
coeff[p1][p] = factorial(p)/factorial(coeffx[p1][p])/factorial(coeffy[p1][p])/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];
|
||||
// jnum = numneigh[i];
|
||||
//select the more efficient algorithm for this particular potential and environment.
|
||||
if (jnum*2>(m+1)*m*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];
|
||||
// double **x = sim->x;
|
||||
int *type = sim->type;
|
||||
double cutmax = pair->cutmax;
|
||||
int res = pair->res;
|
||||
double cutinv2 = 1/cutmax/cutmax;
|
||||
ilist = sim->ilist;
|
||||
// numneigh = sim->numneigh;
|
||||
// firstneigh = sim->firstneigh;
|
||||
int nelements=pair->nelements;
|
||||
i = ilist[ii];
|
||||
itype = pair->map[type[i]];
|
||||
int f = pair->net[itype].dimensions[0];
|
||||
// xtmp = x[i][0];
|
||||
// ytmp = x[i][1];
|
||||
// ztmp = x[i][2];
|
||||
// jlist = firstneigh[i];
|
||||
// jnum = numneigh[i];
|
||||
double expr[jnum][this->k+12];
|
||||
int p = this->k;
|
||||
int countmb=((this->m)*(this->m+1))>>1;
|
||||
// calculate interpolation expr, rinvs and dfc, for each neighbor
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[1] != nelements && atomtypes[1] != jtype && atomtypes[2] != nelements && atomtypes[2] != jtype){
|
||||
expr[jj][0]=0;
|
||||
continue;
|
||||
}
|
||||
// delx = xtmp - x[j][0];
|
||||
// dely = ytmp - x[j][1];
|
||||
// delz = ztmp - x[j][2];
|
||||
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)*this->k];
|
||||
double *p1 = &expcuttable[m1*this->k];
|
||||
double *p2 = &expcuttable[(m1+1)*this->k];
|
||||
double *p3 = &expcuttable[(m1+2)*this->k];
|
||||
for (kk=0;kk<this->k;kk++){
|
||||
expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
|
||||
}
|
||||
double* q = &dfctable[m1-1];
|
||||
double* ri = &rinvsqrttable[m1-1];
|
||||
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 rinvs = ri[1] + 0.5 * r1*(ri[2] - ri[0] + r1*(2.0*ri[0] - 5.0*ri[1] + 4.0*ri[2] - ri[3] + r1*(3.0*(ri[1] - ri[2]) + ri[3] - ri[0])));
|
||||
|
||||
expr[jj][p]=delx*rinvs;
|
||||
expr[jj][p+1]=dely*rinvs;
|
||||
expr[jj][p+2]=delz*rinvs;
|
||||
//Hack to avoid nan when x y or z component of radial vector is exactly 0. Shouldn't affect accuracy.
|
||||
if (expr[jj][p]*expr[jj][p]<0.000000000001){
|
||||
expr[jj][p] = 0.000001;
|
||||
}
|
||||
if (expr[jj][p+1]*expr[jj][p+1]<0.000000000001){
|
||||
expr[jj][p+1] = 0.000001;
|
||||
}
|
||||
if (expr[jj][p+2]*expr[jj][p+2]<0.000000000001){
|
||||
expr[jj][p+2] = 0.000001;
|
||||
}
|
||||
expr[jj][p+3] = -dfc*expr[jj][p];
|
||||
expr[jj][p+4] = rinvs/expr[jj][p];
|
||||
expr[jj][p+5] = rinvs*expr[jj][p];
|
||||
expr[jj][p+6] = -dfc*expr[jj][p+1];
|
||||
expr[jj][p+7] = rinvs/expr[jj][p+1];
|
||||
expr[jj][p+8] = rinvs*expr[jj][p+1];
|
||||
expr[jj][p+9] = -dfc*expr[jj][p+2];
|
||||
expr[jj][p+10] = rinvs/expr[jj][p+2];
|
||||
expr[jj][p+11] = rinvs*expr[jj][p+2];
|
||||
}
|
||||
|
||||
int kb = this->k;
|
||||
int mb = this->m;
|
||||
count = startingneuron;
|
||||
double Bb[mb];
|
||||
double dBbx;
|
||||
double dBby;
|
||||
double dBbz;
|
||||
// double dBbx1[mb];
|
||||
// double dBby1[mb];
|
||||
// double dBbz1[mb];
|
||||
double yprod;
|
||||
for (mcount=0;mcount<countmb;mcount++){
|
||||
int *M = Mf[mcount];
|
||||
int *coeffx = this->coeffx[mcount];
|
||||
int *coeffy = this->coeffy[mcount];
|
||||
int *coeffz = this->coeffz[mcount];
|
||||
int *coeff = this->coeff[mcount];
|
||||
a = mb+1;
|
||||
for (a1=0;a1<mb;a1++){
|
||||
if (Mf[mcount][a1+1]==0){
|
||||
a = a1;
|
||||
break;
|
||||
}
|
||||
}
|
||||
for (n=0;n<kb;n++){
|
||||
for (a1=0;a1<mb;a1++){
|
||||
Bb[a1]=0;
|
||||
// dBbx1[a1] = 0;
|
||||
// dBby1[a1] = 0;
|
||||
// dBbz1[a1] = 0;
|
||||
}
|
||||
ai = n;
|
||||
double y1 = alpha_k[ai]/re;
|
||||
//loop over jtype to get Bb
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (expr[jj][0]==0){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[1] != nelements && atomtypes[1] != jtype){
|
||||
continue;
|
||||
}
|
||||
double yprod = expr[jj][ai];
|
||||
double *y4 = &expr[jj][p];
|
||||
for (a2=0;a2<a;a2++){
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
for (a2=a;a2<mb;a2++){
|
||||
Bb[a2]=Bb[a2]+yprod;
|
||||
// if (i==5 && a2 == 1 && n==0){
|
||||
// char str[128];
|
||||
// sprintf(str,"%f %f %f\n",Bb[a2],yprod,y4[M[a2]]);
|
||||
// std::cout<<str;
|
||||
// }
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
}
|
||||
if (atomtypes[1]!=atomtypes[2]){//Bb!=Bg
|
||||
double Bg[mb];
|
||||
// double dBgx;
|
||||
// double dBgy;
|
||||
// double dBgz;
|
||||
// double dBgx1[mb];
|
||||
// double dBgy1[mb];
|
||||
// double dBgz1[mb];
|
||||
for (a1=0;a1<mb;a1++){
|
||||
Bg[a1]=0;
|
||||
// dBgx1[a1] = 0;
|
||||
// dBgy1[a1] = 0;
|
||||
// dBgz1[a1] = 0;
|
||||
}
|
||||
ai = n;
|
||||
double y1 = alpha_k[ai]/re;
|
||||
//loop over ktype to get Bg
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (expr[jj][0]==0){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[2] != nelements && atomtypes[2] != jtype){
|
||||
continue;
|
||||
}
|
||||
double yprod = expr[jj][ai];
|
||||
double *y4 = &expr[jj][p];
|
||||
for (a2=0;a2<a;a2++){
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
for (a2=a;a2<mb;a2++){
|
||||
Bg[a2]=Bg[a2]+yprod;
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
}
|
||||
double B1;
|
||||
//loop over ktype to get dBg*Bb
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (expr[jj][0]==0){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[2] != nelements && atomtypes[2] != jtype){
|
||||
continue;
|
||||
}
|
||||
double *y3 = &expr[jj][p+3];
|
||||
double *y4 = &expr[jj][p];
|
||||
ai = n;
|
||||
yprod = expr[jj][ai];
|
||||
for (a2=0;a2<a;a2++){
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
ai = n*(mb)+a+count+jj*f;
|
||||
for (a2=a;a2<mb;a2++){
|
||||
B1 = Bb[a2]*coeff[a2]*yprod;
|
||||
dBbx = -B1*(y1*y4[0]+y3[0]-coeffx[a2]*y3[1]+a2*y3[2]);
|
||||
// if (j==1 && i==5 && a2 < 2 && n==0){
|
||||
// char str[128];
|
||||
// sprintf(str,"%d %f %f %f %f %d %f %f %f %d %f %f\n\n",a2,dBbx,y1,y4[0],y3[0],coeffx[a2],y3[1],y3[2],Bb[a2],coeff[a2],yprod,xtmp - x[j][0]);
|
||||
// std::cout<<str;
|
||||
// }
|
||||
dBby = -B1*(y1*y4[1]+y3[3]-coeffy[a2]*y3[4]+a2*y3[5]);
|
||||
dBbz = -B1*(y1*y4[2]+y3[6]-coeffz[a2]*y3[7]+a2*y3[8]);
|
||||
// dBbx1[a2] -= dBbx;
|
||||
// dBby1[a2] -= dBby;
|
||||
// dBbz1[a2] -= dBbz;
|
||||
dfeaturesx[ai] += dBbx;
|
||||
dfeaturesy[ai] += dBby;
|
||||
dfeaturesz[ai] += dBbz;
|
||||
yprod *= y4[M[a2+1]];
|
||||
ai++;
|
||||
}
|
||||
}
|
||||
//loop over jtype to get dBb*Bg
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (expr[jj][0]==0){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[1] != nelements && atomtypes[1] != jtype){
|
||||
continue;
|
||||
}
|
||||
double *y3 = &expr[jj][p+3];
|
||||
double *y4 = &expr[jj][p];
|
||||
ai = n;
|
||||
yprod = expr[jj][ai];
|
||||
for (a2=0;a2<a;a2++){
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
ai = n*(mb)+a+count+jj*f;
|
||||
for (a2=a;a2<mb;a2++){
|
||||
B1 = Bg[a2]*coeff[a2]*yprod;
|
||||
dBbx = -B1*(y1*y4[0]+y3[0]-coeffx[a2]*y3[1]+a2*y3[2]);
|
||||
dBby = -B1*(y1*y4[1]+y3[3]-coeffy[a2]*y3[4]+a2*y3[5]);
|
||||
dBbz = -B1*(y1*y4[2]+y3[6]-coeffz[a2]*y3[7]+a2*y3[8]);
|
||||
// dBbx1[a2] -= dBbx;
|
||||
// dBby1[a2] -= dBby;
|
||||
// dBbz1[a2] -= dBbz;
|
||||
dfeaturesx[ai] += dBbx;
|
||||
dfeaturesy[ai] += dBby;
|
||||
dfeaturesz[ai] += dBbz;
|
||||
yprod *= y4[M[a2+1]];
|
||||
ai++;
|
||||
}
|
||||
}
|
||||
//central atom derivative
|
||||
for (a2=a;a2<mb;a2++){
|
||||
ai = n*(mb)+a2+count+jnum*f;
|
||||
features[ai-jnum*f] += Bb[a2]*Bg[a2]*coeff[a2];
|
||||
// dfeaturesx[ai] += dBbx1[a2];
|
||||
// dfeaturesy[ai] += dBby1[a2];
|
||||
// dfeaturesz[ai] += dBbz1[a2];
|
||||
}
|
||||
}
|
||||
else{//Bb=Bg
|
||||
double B1;
|
||||
//loop over jtype to get 2*Bb*dBb
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (expr[jj][0]==0){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[1] != nelements && atomtypes[1] != jtype){
|
||||
continue;
|
||||
}
|
||||
double *y3 = &expr[jj][p+3];
|
||||
double *y4 = &expr[jj][p];
|
||||
ai = n;
|
||||
yprod = expr[jj][ai];
|
||||
for (a2=0;a2<a;a2++){
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
ai = n*(mb)+a+count+jj*f;
|
||||
for (a2=a;a2<mb;a2++){
|
||||
B1 = 2*Bb[a2]*coeff[a2]*yprod;
|
||||
dBbx = -B1*(y1*y4[0]+y3[0]-coeffx[a2]*y3[1]+a2*y3[2]);
|
||||
dBby = -B1*(y1*y4[1]+y3[3]-coeffy[a2]*y3[4]+a2*y3[5]);
|
||||
dBbz = -B1*(y1*y4[2]+y3[6]-coeffz[a2]*y3[7]+a2*y3[8]);
|
||||
// dBbx1[a2] -= dBbx;
|
||||
// dBby1[a2] -= dBby;
|
||||
// dBbz1[a2] -= dBbz;
|
||||
dfeaturesx[ai] += dBbx;
|
||||
dfeaturesy[ai] += dBby;
|
||||
dfeaturesz[ai] += dBbz;
|
||||
yprod *= y4[M[a2+1]];
|
||||
ai++;
|
||||
}
|
||||
}
|
||||
//central atom derivative
|
||||
for (a2=a;a2<mb;a2++){
|
||||
ai = n*(mb)+a2+count+jnum*f;
|
||||
features[ai-jnum*f] += Bb[a2]*Bb[a2]*coeff[a2];
|
||||
// dfeaturesx[ai] += dBbx1[a2];
|
||||
// dfeaturesy[ai] += dBby1[a2];
|
||||
// dfeaturesz[ai] += dBbz1[a2];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (expr[jj][0]==0){continue;}
|
||||
count = startingneuron;
|
||||
for (n=0;n<kb;n++){
|
||||
for (m=0;m<mb;m++){
|
||||
dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
|
||||
dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
|
||||
dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
|
||||
count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//Called by do3bodyfeatureset. Algorithm for low neighbor numbers and large series of bond angle powers
|
||||
void Fingerprint_bond::do3bodyfeatureset_doubleneighborloop(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,ktype,kk,m,n;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
//int itypes = atomtypes[0];
|
||||
int jtypes = atomtypes[1];
|
||||
int ktypes = atomtypes[2];
|
||||
int count=0;
|
||||
PairRANN::Simulation *sim = &pair->sims[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;
|
||||
// numneigh = sim->numneigh;
|
||||
// firstneigh = sim->firstneigh;
|
||||
i = ilist[ii];
|
||||
itype = pair->map[type[i]];
|
||||
int f = pair->net[itype].dimensions[0];
|
||||
// xtmp = x[i][0];
|
||||
// ytmp = x[i][1];
|
||||
// ztmp = x[i][2];
|
||||
// jlist = firstneigh[i];
|
||||
// jnum = numneigh[i];
|
||||
double expr[jnum][this->k];
|
||||
double y[jnum][3];
|
||||
double ri[jnum];
|
||||
double dfc[jnum];
|
||||
int kb = this->k;
|
||||
int mb = this->m;
|
||||
double c41[this->k];
|
||||
double c51[this->k];
|
||||
double c61[this->k];
|
||||
double ct[this->k];
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype){
|
||||
expr[jj][0]=0;
|
||||
continue;
|
||||
}
|
||||
// delx = xtmp - x[j][0];
|
||||
// dely = ytmp - x[j][1];
|
||||
// delz = ztmp - x[j][2];
|
||||
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("Neighbor list is invalid.");//usually results from nan somewhere.
|
||||
r1 = r1-trunc(r1);
|
||||
double *p0 = &expcuttable[(m1-1)*this->k];
|
||||
double *p1 = &expcuttable[m1*this->k];
|
||||
double *p2 = &expcuttable[(m1+1)*this->k];
|
||||
double *p3 = &expcuttable[(m1+2)*this->k];
|
||||
for (kk=0;kk<this->k;kk++){
|
||||
expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
|
||||
}
|
||||
double* q = &dfctable[m1-1];
|
||||
double* r2 = &rinvsqrttable[m1-1];
|
||||
dfc[jj] = 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])));
|
||||
ri[jj] = r2[1] + 0.5 * r1*(r2[2] - r2[0] + r1*(2.0*r2[0] - 5.0*r2[1] + 4.0*r2[2] - r2[3] + r1*(3.0*(r2[1] - r2[2]) + r2[3] - r2[0])));
|
||||
y[jj][0]=delx*ri[jj];
|
||||
y[jj][1]=dely*ri[jj];
|
||||
y[jj][2]=delz*ri[jj];
|
||||
}
|
||||
// if (i==5){
|
||||
// for (jj=0;jj<jnum;jj++){
|
||||
// j = jlist[jj];
|
||||
// if (j==45){
|
||||
// sprintf(str,"%d %f %f %f\n",jj,dfeaturesx[jj*f+0],dfeaturesy[jj*f+0],dfeaturesz[jj*f+0]);
|
||||
// std::cout<<str;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
if (expr[jj][0]==0)continue;
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (jtypes != nelements && jtypes != jtype){
|
||||
continue;
|
||||
}
|
||||
for (n = 0;n<this->k;n++){
|
||||
ct[n] = 2*alpha_k[n]/re;
|
||||
c41[n]=(-ct[n]+2*dfc[jj])*y[jj][0];
|
||||
c51[n]=(-ct[n]+2*dfc[jj])*y[jj][1];
|
||||
c61[n]= (-ct[n]+2*dfc[jj])*y[jj][2];
|
||||
}
|
||||
if (jtypes==ktypes){
|
||||
for (kk=jj+1;kk< jnum; kk++){
|
||||
if (expr[kk][0]==0)continue;
|
||||
// int k1 = jlist[kk];
|
||||
// k1 &= NEIGHMASK;
|
||||
// ktype = pair->map[type[k1]];
|
||||
ktype = tn[kk];
|
||||
if (ktypes != nelements && ktypes != ktype){
|
||||
continue;
|
||||
}
|
||||
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]);
|
||||
// double c1 = 2*ri[jj]*y[kk][0]*(1-y[jj][0]*y[jj][0]);
|
||||
// double c2 = 2*ri[jj]*y[kk][1]*(1-y[jj][1]*y[jj][1]);
|
||||
// double c3 = 2*ri[jj]*y[kk][2]*(1-y[jj][2]*y[jj][2]);
|
||||
// double c10 = 2*ri[kk]*y[jj][0]*(1-y[kk][0]*y[kk][0]);
|
||||
// double c11 = 2*ri[kk]*y[jj][1]*(1-y[kk][1]*y[kk][1]);
|
||||
// double c12 = 2*ri[kk]*y[jj][2]*(1-y[kk][2]*y[kk][2]);
|
||||
for (n=0;n<kb;n++){
|
||||
double dot1=expr[jj][n]*expr[kk][n];
|
||||
double c4 = c41[n];
|
||||
double c5 = c51[n];
|
||||
double c6 = c61[n];
|
||||
double ct2 = -ct[n]+2*dfc[kk];
|
||||
double c42 = ct2*y[kk][0];
|
||||
double c52 = ct2*y[kk][1];
|
||||
double c62 = ct2*y[kk][2];
|
||||
//m=0
|
||||
features[count]+=2*dot1;
|
||||
dfeaturesx[jj*f+count]+=dot1*c4;
|
||||
dfeaturesy[jj*f+count]+=dot1*c5;
|
||||
dfeaturesz[jj*f+count]+=dot1*c6;
|
||||
dfeaturesx[kk*f+count]+=dot1*c42;
|
||||
dfeaturesy[kk*f+count]+=dot1*c52;
|
||||
dfeaturesz[kk*f+count]+=dot1*c62;
|
||||
c4*=dot;
|
||||
c5*=dot;
|
||||
c6*=dot;
|
||||
c42*=dot;
|
||||
c52*=dot;
|
||||
c62*=dot;
|
||||
count++;
|
||||
for (m=1;m<mb;m++){
|
||||
double c7 = dot1*(m*c1+c4);
|
||||
double c8 = dot1*(m*c2+c5);
|
||||
double c9 = dot1*(m*c3+c6);
|
||||
dfeaturesx[jj*f+count]+=c7;
|
||||
dfeaturesy[jj*f+count]+=c8;
|
||||
dfeaturesz[jj*f+count]+=c9;
|
||||
dfeaturesx[kk*f+count]+=dot1*(m*c10+c42);
|
||||
dfeaturesy[kk*f+count]+=dot1*(m*c11+c52);
|
||||
dfeaturesz[kk*f+count]+=dot1*(m*c12+c62);
|
||||
dot1*=dot;
|
||||
features[count++]+=2*dot1;
|
||||
}
|
||||
}
|
||||
}
|
||||
kk=jj;
|
||||
if (ktypes == nelements || ktypes == jtype){
|
||||
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 c1 = 2*ri[jj]*y[kk][0]*(1-y[jj][0]*y[jj][0]);
|
||||
// double c2 = 2*ri[jj]*y[kk][1]*(1-y[jj][1]*y[jj][1]);
|
||||
// double c3 = 2*ri[jj]*y[kk][2]*(1-y[jj][2]*y[jj][2]);
|
||||
for (n=0;n<kb;n++){
|
||||
double dot1=expr[jj][n]*expr[kk][n];
|
||||
double c4 = c41[n];
|
||||
double c5 = c51[n];
|
||||
double c6 = c61[n];
|
||||
//m=0
|
||||
features[count]+=dot1;
|
||||
dfeaturesx[jj*f+count]+=dot1*c4;
|
||||
dfeaturesy[jj*f+count]+=dot1*c5;
|
||||
dfeaturesz[jj*f+count]+=dot1*c6;
|
||||
c4*=dot;
|
||||
c5*=dot;
|
||||
c6*=dot;
|
||||
count++;
|
||||
for (m=1;m<mb;m++){
|
||||
double c7 = dot1*(m*c1+c4);
|
||||
double c8 = dot1*(m*c2+c5);
|
||||
double c9 = dot1*(m*c3+c6);
|
||||
dfeaturesx[jj*f+count]+=c7;
|
||||
dfeaturesy[jj*f+count]+=c8;
|
||||
dfeaturesz[jj*f+count]+=c9;
|
||||
dot1*=dot;
|
||||
features[count++]+=dot1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
for (kk=0;kk<jnum; kk++){
|
||||
if (expr[kk][0]==0)continue;
|
||||
// int k1 = jlist[kk];
|
||||
// k1 &= NEIGHMASK;
|
||||
// ktype = pair->map[type[k1]];
|
||||
ktype = tn[kk];
|
||||
if (ktypes != nelements && ktypes != ktype){
|
||||
continue;
|
||||
}
|
||||
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]);
|
||||
// double c1 = 2*ri[jj]*y[kk][0]*(1-y[jj][0]*y[jj][0]);
|
||||
// double c2 = 2*ri[jj]*y[kk][1]*(1-y[jj][1]*y[jj][1]);
|
||||
// double c3 = 2*ri[jj]*y[kk][2]*(1-y[jj][2]*y[jj][2]);
|
||||
// double c10 = 2*ri[kk]*y[jj][0]*(1-y[kk][0]*y[kk][0]);
|
||||
// double c11 = 2*ri[kk]*y[jj][1]*(1-y[kk][1]*y[kk][1]);
|
||||
// double c12 = 2*ri[kk]*y[jj][2]*(1-y[kk][2]*y[kk][2]);
|
||||
for (n=0;n<kb;n++){
|
||||
double dot1=expr[jj][n]*expr[kk][n];
|
||||
double c4 = c41[n]/2;
|
||||
double c5 = c51[n]/2;
|
||||
double c6 = c61[n]/2;
|
||||
double ct2 = -ct[n]/2+dfc[kk];
|
||||
double c42 = ct2*y[kk][0];
|
||||
double c52 = ct2*y[kk][1];
|
||||
double c62 = ct2*y[kk][2];
|
||||
//m=0
|
||||
features[count]+=dot1;
|
||||
dfeaturesx[jj*f+count]+=dot1*c4;
|
||||
dfeaturesy[jj*f+count]+=dot1*c5;
|
||||
dfeaturesz[jj*f+count]+=dot1*c6;
|
||||
dfeaturesx[kk*f+count]+=dot1*c42;
|
||||
dfeaturesy[kk*f+count]+=dot1*c52;
|
||||
dfeaturesz[kk*f+count]+=dot1*c62;
|
||||
c4*=dot;
|
||||
c5*=dot;
|
||||
c6*=dot;
|
||||
c42*=dot;
|
||||
c52*=dot;
|
||||
c62*=dot;
|
||||
count++;
|
||||
for (m=1;m<mb;m++){
|
||||
double c7 = dot1*(m*c1+c4);
|
||||
double c8 = dot1*(m*c2+c5);
|
||||
double c9 = dot1*(m*c3+c6);
|
||||
dfeaturesx[jj*f+count]+=c7;
|
||||
dfeaturesy[jj*f+count]+=c8;
|
||||
dfeaturesz[jj*f+count]+=c9;
|
||||
dfeaturesx[kk*f+count]+=dot1*(m*c10+c42);
|
||||
dfeaturesy[kk*f+count]+=dot1*(m*c11+c52);
|
||||
dfeaturesz[kk*f+count]+=dot1*(m*c12+c62);
|
||||
dot1*=dot;
|
||||
features[count++]+=dot1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (expr[jj][0]==0){continue;}
|
||||
count = startingneuron;
|
||||
for (n=0;n<kb;n++){
|
||||
for (m=0;m<mb;m++){
|
||||
dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
|
||||
dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
|
||||
dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
|
||||
count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int Fingerprint_bond::factorial(int n) {
|
||||
if ((n==0)||(n==1))
|
||||
return 1;
|
||||
else
|
||||
return n*factorial(n-1);
|
||||
}
|
||||
67
src/fingerprint_bond.h
Normal file
67
src/fingerprint_bond.h
Normal file
@ -0,0 +1,67 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#ifdef FINGERPRINT_CLASS
|
||||
|
||||
FingerprintStyle(bond,Fingerprint_bond)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef FINGERPRINT_BOND_H_
|
||||
#define FINGERPRINT_BOND_H_
|
||||
|
||||
#include "fingerprint.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Fingerprint_bond : public Fingerprint {
|
||||
public:
|
||||
Fingerprint_bond(PairRANN *);
|
||||
~Fingerprint_bond();
|
||||
bool parse_values(char*,char*);
|
||||
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 *);
|
||||
int factorial(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 k;
|
||||
int m;
|
||||
int **Mf;
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
#endif /* FINGERPRINT_BOND_H_ */
|
||||
946
src/fingerprint_bondscreened.cpp
Normal file
946
src/fingerprint_bondscreened.cpp
Normal file
@ -0,0 +1,946 @@
|
||||
/* -*- 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@cavs.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
|
||||
#include "fingerprint_bondscreened.h"
|
||||
#include "fingerprint.h"
|
||||
#include <math.h>
|
||||
#include <cmath>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <iostream>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
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;
|
||||
k = 0;
|
||||
m = 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<(m*(m+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(char * constant, char * line1){
|
||||
char **words=new char *[MAXLINE];
|
||||
int nwords,l;
|
||||
nwords=0;
|
||||
words[nwords++] = strtok(line1,": ,\t\n");
|
||||
while ((words[nwords++] = strtok(NULL,": ,\t\n"))) continue;
|
||||
nwords -= 1;
|
||||
if (strcmp(constant,"re")==0){
|
||||
re = strtod(words[0],NULL);
|
||||
}
|
||||
else if (strcmp(constant,"rc")==0){
|
||||
rc = strtod(words[0],NULL);
|
||||
}
|
||||
else if (strcmp(constant,"alphak")==0){
|
||||
delete [] alpha_k;
|
||||
alpha_k = new double [nwords];
|
||||
for (l=0;l<nwords;l++){
|
||||
alpha_k[l]=strtod(words[l],NULL);
|
||||
}
|
||||
}
|
||||
else if (strcmp(constant,"dr")==0){
|
||||
dr = strtod(words[0],NULL);
|
||||
}
|
||||
else if (strcmp(constant,"k")==0){
|
||||
k = strtol(words[0],NULL,10);
|
||||
}
|
||||
else if (strcmp(constant,"m")==0){
|
||||
m = strtol(words[0],NULL,10);
|
||||
}
|
||||
else pair->errorf("Undefined value for bond power");
|
||||
delete [] words;
|
||||
if (re!=0.0 && rc!=0.0 && alpha_k[0]!=-1 && dr!=0.0 && m!=0 && k!=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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:alphak:\n",style,id);
|
||||
for (i=0;i<k;i++){
|
||||
fprintf(fid,"%f ",alpha_k[i]);
|
||||
}
|
||||
fprintf(fid,"\n");
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:k:\n",style,id);
|
||||
fprintf(fid,"%d\n",k);
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:m:\n",style,id);
|
||||
fprintf(fid,"%d\n",m);
|
||||
}
|
||||
|
||||
void Fingerprint_bondscreened::init(int *i,int id){
|
||||
for (int j=0;j<n_body_type;j++){atomtypes[j] = i[j];}
|
||||
re = 0;
|
||||
rc = 0;
|
||||
m = 0;
|
||||
this->k = 0;
|
||||
alpha_k = new double [1];
|
||||
alpha_k[0]=-1;
|
||||
empty = false;
|
||||
this->id = id;
|
||||
}
|
||||
|
||||
//number of neurons defined by this fingerprint
|
||||
int Fingerprint_bondscreened::get_length(){
|
||||
return m*k;
|
||||
}
|
||||
|
||||
void Fingerprint_bondscreened::allocate(){
|
||||
generate_exp_cut_table();
|
||||
generate_coefficients();
|
||||
generate_rinvssqrttable();
|
||||
}
|
||||
|
||||
//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop and do3bodyfeatureset_doubleneighborloop.
|
||||
void Fingerprint_bondscreened::generate_exp_cut_table(){
|
||||
int m,n;
|
||||
double r1;
|
||||
int buf = 5;
|
||||
int res = pair->res;
|
||||
double cutmax = pair->cutmax;
|
||||
expcuttable = new double [(res+buf)*(this->k)];
|
||||
dfctable = new double [res+buf];
|
||||
for (m=0;m<(res+buf);m++){
|
||||
r1 = cutmax*cutmax*(double)(m)/(double)(res);
|
||||
for (n=0;n<(this->k);n++){
|
||||
expcuttable[n+m*(this->k)] = 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 = this->m;
|
||||
mc=(mb*(mb+1))>>1;
|
||||
coeff = new int *[mc];
|
||||
coeffx = new int *[mc];
|
||||
coeffy = new int *[mc];
|
||||
coeffz = new int *[mc];
|
||||
for (p=0;p<mc;p++){
|
||||
coeff[p]=new int [mb];
|
||||
coeffx[p]=new int [mb];
|
||||
coeffy[p]=new int [mb];
|
||||
coeffz[p]=new int [mb];
|
||||
}
|
||||
Mf = new int*[mc];
|
||||
int *M = new int[this->m+1];
|
||||
for (p=0;p<this->m+1;p++){
|
||||
M[p]=0;
|
||||
}
|
||||
for (p1=0;p1<mc;p1++){
|
||||
Mf[p1] = new int[m+1];
|
||||
for (p=0;p<m+1;p++){
|
||||
Mf[p1][p]=0;
|
||||
}
|
||||
}
|
||||
M[0] = 2;
|
||||
Mf[0][0] = 2;
|
||||
n = 1;
|
||||
int m1 = 1;
|
||||
bool go = true;
|
||||
bool broke = false;
|
||||
while (go){
|
||||
broke = false;
|
||||
for (i1=0;i1<m-1;i1++){
|
||||
if (M[i1+1] == 0){
|
||||
M[i1+1]=M[i1+1]+1;
|
||||
for (p1=0;p1<m+1;p1++){
|
||||
Mf[n][p1] = M[p1];
|
||||
}
|
||||
n = n+1;
|
||||
broke = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (m1<m && !broke){
|
||||
M[m1]=M[m1]+1;
|
||||
for (p1=m1+1;p1<m+1;p1++){
|
||||
M[p1]=0;
|
||||
}
|
||||
for (p1=0;p1<m+1;p1++){
|
||||
Mf[n][p1]=M[p1];
|
||||
}
|
||||
n=n+1;
|
||||
broke = true;
|
||||
m1 = m1+1;
|
||||
}
|
||||
if (!broke){
|
||||
go = false;
|
||||
}
|
||||
}
|
||||
for (p=0;p<mb;p++){
|
||||
for (p1=0;p1<mc;p1++){
|
||||
if (p==0){
|
||||
coeffx[p1][p]=0;
|
||||
coeffy[p1][p]=0;
|
||||
coeffz[p1][p]=0;
|
||||
}
|
||||
else{
|
||||
coeffx[p1][p]=coeffx[p1][p-1];
|
||||
if (Mf[p1][p]==0){
|
||||
coeffx[p1][p]++;
|
||||
}
|
||||
coeffy[p1][p]=coeffy[p1][p-1];
|
||||
if (Mf[p1][p]==1){
|
||||
coeffy[p1][p]++;
|
||||
}
|
||||
coeffz[p1][p]=coeffz[p1][p-1];
|
||||
if (Mf[p1][p]==2){
|
||||
coeffz[p1][p]++;
|
||||
}
|
||||
}
|
||||
coeff[p1][p] = factorial(p)/factorial(coeffx[p1][p])/factorial(coeffy[p1][p])/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;
|
||||
// int *ilist,*numneigh;
|
||||
// PairRANN::Simulation *sim = &pair->sims[sid];
|
||||
// ilist = sim->ilist;
|
||||
// numneigh = sim->numneigh;
|
||||
// i = ilist[ii];
|
||||
// jnum = numneigh[i];
|
||||
//select the more efficient algorithm for this particular potential and environment.
|
||||
if (jnum*2>(m+1)*m*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];
|
||||
// double **x = sim->x;
|
||||
//double **f = atom->f;
|
||||
int *type = sim->type;
|
||||
double cutmax = pair->cutmax;
|
||||
int res = pair->res;
|
||||
double cutinv2 = 1/cutmax/cutmax;
|
||||
ilist = sim->ilist;
|
||||
// numneigh = sim->numneigh;
|
||||
// firstneigh = sim->firstneigh;
|
||||
int nelements=pair->nelements;
|
||||
i = ilist[ii];
|
||||
itype = pair->map[type[i]];
|
||||
int f = pair->net[itype].dimensions[0];
|
||||
// xtmp = x[i][0];
|
||||
// ytmp = x[i][1];
|
||||
// ztmp = x[i][2];
|
||||
// jlist = firstneigh[i];
|
||||
// jnum = numneigh[i];
|
||||
double expr[jnum][this->k+12];
|
||||
int p = this->k;
|
||||
int countmb=((this->m)*(this->m+1))>>1;
|
||||
// calculate interpolation expr, rinvs and dfc, for each neighbor
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
if (Bij[jj]==false){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[1] != nelements && atomtypes[1] != jtype && atomtypes[2] != nelements && atomtypes[2] != jtype){
|
||||
expr[jj][0]=0;
|
||||
continue;
|
||||
}
|
||||
// delx = xtmp - x[j][0];
|
||||
// dely = ytmp - x[j][1];
|
||||
// delz = ztmp - x[j][2];
|
||||
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)*this->k];
|
||||
double *p1 = &expcuttable[m1*this->k];
|
||||
double *p2 = &expcuttable[(m1+1)*this->k];
|
||||
double *p3 = &expcuttable[(m1+2)*this->k];
|
||||
for (kk=0;kk<this->k;kk++){
|
||||
expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
|
||||
expr[jj][kk] *= Sik[jj];
|
||||
}
|
||||
double* q = &dfctable[m1-1];
|
||||
double* ri = &rinvsqrttable[m1-1];
|
||||
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 rinvs = ri[1] + 0.5 * r1*(ri[2] - ri[0] + r1*(2.0*ri[0] - 5.0*ri[1] + 4.0*ri[2] - ri[3] + r1*(3.0*(ri[1] - ri[2]) + ri[3] - ri[0])));
|
||||
|
||||
expr[jj][p]=delx*rinvs;
|
||||
expr[jj][p+1]=dely*rinvs;
|
||||
expr[jj][p+2]=delz*rinvs;
|
||||
//Hack to avoid nan when x y or z component of radial vector is exactly 0. Shouldn't affect accuracy.
|
||||
if (expr[jj][p]*expr[jj][p]<0.000000000001){
|
||||
expr[jj][p] = 0.000001;
|
||||
}
|
||||
if (expr[jj][p+1]*expr[jj][p+1]<0.000000000001){
|
||||
expr[jj][p+1] = 0.000001;
|
||||
}
|
||||
if (expr[jj][p+2]*expr[jj][p+2]<0.000000000001){
|
||||
expr[jj][p+2] = 0.000001;
|
||||
}
|
||||
expr[jj][p+3] = -dfc*expr[jj][p]-dSikx[jj];
|
||||
expr[jj][p+4] = rinvs/expr[jj][p];
|
||||
expr[jj][p+5] = rinvs*expr[jj][p];
|
||||
expr[jj][p+6] = -dfc*expr[jj][p+1]-dSiky[jj];
|
||||
expr[jj][p+7] = rinvs/expr[jj][p+1];
|
||||
expr[jj][p+8] = rinvs*expr[jj][p+1];
|
||||
expr[jj][p+9] = -dfc*expr[jj][p+2]-dSikz[jj];
|
||||
expr[jj][p+10] = rinvs/expr[jj][p+2];
|
||||
expr[jj][p+11] = rinvs*expr[jj][p+2];
|
||||
}
|
||||
|
||||
int kb = this->k;
|
||||
int mb = this->m;
|
||||
// int ind = kb+mb*kb;
|
||||
count = startingneuron;
|
||||
double Bb[mb];
|
||||
double dBbx;
|
||||
double dBby;
|
||||
double dBbz;
|
||||
double dBbx1[mb];
|
||||
double dBby1[mb];
|
||||
double dBbz1[mb];
|
||||
double yprod;
|
||||
for (mcount=0;mcount<countmb;mcount++){
|
||||
int *M = Mf[mcount];
|
||||
int *coeffx = this->coeffx[mcount];
|
||||
int *coeffy = this->coeffy[mcount];
|
||||
int *coeffz = this->coeffz[mcount];
|
||||
int *coeff = this->coeff[mcount];
|
||||
a = mb+1;
|
||||
for (a1=0;a1<mb;a1++){
|
||||
if (Mf[mcount][a1+1]==0){
|
||||
a = a1;
|
||||
break;
|
||||
}
|
||||
}
|
||||
for (n=0;n<kb;n++){
|
||||
for (a1=0;a1<mb;a1++){
|
||||
Bb[a1]=0;
|
||||
dBbx1[a1] = 0;
|
||||
dBby1[a1] = 0;
|
||||
dBbz1[a1] = 0;
|
||||
}
|
||||
ai = n;
|
||||
double y1 = alpha_k[ai]/re;
|
||||
//loop over jtype to get Bb
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (Bij[jj]==false){continue;}
|
||||
if (expr[jj][0]==0){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[1] != nelements && atomtypes[1] != jtype){
|
||||
continue;
|
||||
}
|
||||
double yprod = expr[jj][ai];
|
||||
double *y4 = &expr[jj][p];
|
||||
for (a2=0;a2<a;a2++){
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
for (a2=a;a2<mb;a2++){
|
||||
Bb[a2]=Bb[a2]+yprod;
|
||||
// if (i==5 && a2 == 1 && n==0){
|
||||
// char str[128];
|
||||
// sprintf(str,"%f %f %f\n",Bb[a2],yprod,y4[M[a2]]);
|
||||
// std::cout<<str;
|
||||
// }
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
}
|
||||
if (atomtypes[1]!=atomtypes[2]){//Bb!=Bg
|
||||
double Bg[mb];
|
||||
// double dBgx;
|
||||
// double dBgy;
|
||||
// double dBgz;
|
||||
// double dBgx1[mb];
|
||||
// double dBgy1[mb];
|
||||
// double dBgz1[mb];
|
||||
for (a1=0;a1<mb;a1++){
|
||||
Bg[a1]=0;
|
||||
// dBgx1[a1] = 0;
|
||||
// dBgy1[a1] = 0;
|
||||
// dBgz1[a1] = 0;
|
||||
}
|
||||
ai = n;
|
||||
double y1 = alpha_k[ai]/re;
|
||||
//loop over ktype to get Bg
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (Bij[jj]==false){continue;}
|
||||
if (expr[jj][0]==0){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[2] != nelements && atomtypes[2] != jtype){
|
||||
continue;
|
||||
}
|
||||
double yprod = expr[jj][ai];
|
||||
double *y4 = &expr[jj][p];
|
||||
for (a2=0;a2<a;a2++){
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
for (a2=a;a2<mb;a2++){
|
||||
Bg[a2]=Bg[a2]+yprod;
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
}
|
||||
double B1;
|
||||
//loop over ktype to get dBg*Bb
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (Bij[jj]==false){continue;}
|
||||
if (expr[jj][0]==0){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[2] != nelements && atomtypes[2] != jtype){
|
||||
continue;
|
||||
}
|
||||
double *y3 = &expr[jj][p+3];
|
||||
double *y4 = &expr[jj][p];
|
||||
ai = n;
|
||||
yprod = expr[jj][ai];
|
||||
for (a2=0;a2<a;a2++){
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
ai = n*(mb)+a+count+jj*f;
|
||||
for (a2=a;a2<mb;a2++){
|
||||
B1 = Bb[a2]*coeff[a2]*yprod;
|
||||
dBbx = -B1*(y1*y4[0]+y3[0]-coeffx[a2]*y3[1]+a2*y3[2]);
|
||||
dBby = -B1*(y1*y4[1]+y3[3]-coeffy[a2]*y3[4]+a2*y3[5]);
|
||||
dBbz = -B1*(y1*y4[2]+y3[6]-coeffz[a2]*y3[7]+a2*y3[8]);
|
||||
dBbx1[a2] -= dBbx;
|
||||
dBby1[a2] -= dBby;
|
||||
dBbz1[a2] -= dBbz;
|
||||
dfeaturesx[ai] += dBbx;
|
||||
dfeaturesy[ai] += dBby;
|
||||
dfeaturesz[ai] += dBbz;
|
||||
yprod *= y4[M[a2+1]];
|
||||
ai++;
|
||||
for (kk=0;kk<jnum;kk++){
|
||||
if (Bij[kk]==false){continue;}
|
||||
dfeaturesx[n*mb+a2+count+kk*f]+=B1*dSijkx[jj*jnum+kk];
|
||||
dfeaturesy[n*mb+a2+count+kk*f]+=B1*dSijky[jj*jnum+kk];
|
||||
dfeaturesz[n*mb+a2+count+kk*f]+=B1*dSijkz[jj*jnum+kk];
|
||||
}
|
||||
}
|
||||
}
|
||||
//loop over jtype to get dBb*Bg
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (Bij[jj]==false){continue;}
|
||||
if (expr[jj][0]==0){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[1] != nelements && atomtypes[1] != jtype){
|
||||
continue;
|
||||
}
|
||||
double *y3 = &expr[jj][p+3];
|
||||
double *y4 = &expr[jj][p];
|
||||
ai = n;
|
||||
yprod = expr[jj][ai];
|
||||
for (a2=0;a2<a;a2++){
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
ai = n*(mb)+a+count+jj*f;
|
||||
for (a2=a;a2<mb;a2++){
|
||||
B1 = Bg[a2]*coeff[a2]*yprod;
|
||||
dBbx = -B1*(y1*y4[0]+y3[0]-coeffx[a2]*y3[1]+a2*y3[2]);
|
||||
dBby = -B1*(y1*y4[1]+y3[3]-coeffy[a2]*y3[4]+a2*y3[5]);
|
||||
dBbz = -B1*(y1*y4[2]+y3[6]-coeffz[a2]*y3[7]+a2*y3[8]);
|
||||
dBbx1[a2] -= dBbx;
|
||||
dBby1[a2] -= dBby;
|
||||
dBbz1[a2] -= dBbz;
|
||||
dfeaturesx[ai] += dBbx;
|
||||
dfeaturesy[ai] += dBby;
|
||||
dfeaturesz[ai] += dBbz;
|
||||
yprod *= y4[M[a2+1]];
|
||||
ai++;
|
||||
for (kk=0;kk<jnum;kk++){
|
||||
if (Bij[kk]==false){continue;}
|
||||
dfeaturesx[n*mb+a2+count+kk*f]+=B1*dSijkx[jj*jnum+kk];
|
||||
dfeaturesy[n*mb+a2+count+kk*f]+=B1*dSijky[jj*jnum+kk];
|
||||
dfeaturesz[n*mb+a2+count+kk*f]+=B1*dSijkz[jj*jnum+kk];
|
||||
}
|
||||
}
|
||||
}
|
||||
//central atom derivative
|
||||
for (a2=a;a2<mb;a2++){
|
||||
ai = n*(mb)+a2+count+jnum*f;
|
||||
features[ai-jnum*f] += Bb[a2]*Bg[a2]*coeff[a2];
|
||||
// dfeaturesx[ai] += dBbx1[a2];
|
||||
// dfeaturesy[ai] += dBby1[a2];
|
||||
// dfeaturesz[ai] += dBbz1[a2];
|
||||
}
|
||||
}
|
||||
else{//Bb=Bg
|
||||
double B1;
|
||||
//loop over jtype to get 2*Bb*dBb
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (Bij[jj]==false){continue;}
|
||||
if (expr[jj][0]==0){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (atomtypes[1] != nelements && atomtypes[1] != jtype){
|
||||
continue;
|
||||
}
|
||||
double *y3 = &expr[jj][p+3];
|
||||
double *y4 = &expr[jj][p];
|
||||
ai = n;
|
||||
yprod = expr[jj][ai];
|
||||
for (a2=0;a2<a;a2++){
|
||||
yprod *= y4[M[a2+1]];
|
||||
}
|
||||
ai = n*(mb)+a+count+jj*f;
|
||||
for (a2=a;a2<mb;a2++){
|
||||
B1 = 2*Bb[a2]*coeff[a2]*yprod;
|
||||
dBbx = -B1*(y1*y4[0]+y3[0]-coeffx[a2]*y3[1]+a2*y3[2]);
|
||||
dBby = -B1*(y1*y4[1]+y3[3]-coeffy[a2]*y3[4]+a2*y3[5]);
|
||||
dBbz = -B1*(y1*y4[2]+y3[6]-coeffz[a2]*y3[7]+a2*y3[8]);
|
||||
dBbx1[a2] -= dBbx;
|
||||
dBby1[a2] -= dBby;
|
||||
dBbz1[a2] -= dBbz;
|
||||
dfeaturesx[ai] += dBbx;
|
||||
dfeaturesy[ai] += dBby;
|
||||
dfeaturesz[ai] += dBbz;
|
||||
yprod *= y4[M[a2+1]];
|
||||
ai++;
|
||||
for (kk=0;kk<jnum;kk++){
|
||||
if (Bij[kk]==false){continue;}
|
||||
dfeaturesx[n*mb+a2+count+kk*f]+=B1*dSijkx[jj*jnum+kk];
|
||||
dfeaturesy[n*mb+a2+count+kk*f]+=B1*dSijky[jj*jnum+kk];
|
||||
dfeaturesz[n*mb+a2+count+kk*f]+=B1*dSijkz[jj*jnum+kk];
|
||||
}
|
||||
}
|
||||
}
|
||||
//central atom derivative
|
||||
for (a2=a;a2<mb;a2++){
|
||||
ai = n*(mb)+a2+count+jnum*f;
|
||||
features[ai-jnum*f] += Bb[a2]*Bb[a2]*coeff[a2];
|
||||
// dfeaturesx[ai] += dBbx1[a2];
|
||||
// dfeaturesy[ai] += dBby1[a2];
|
||||
// dfeaturesz[ai] += dBbz1[a2];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (Bij[jj]==false){continue;}
|
||||
if (expr[jj][0]==0){continue;}
|
||||
count = startingneuron;
|
||||
for (n=0;n<kb;n++){
|
||||
for (m=0;m<mb;m++){
|
||||
dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
|
||||
dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
|
||||
dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
|
||||
count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//Called by do3bodyfeatureset. Algorithm for low neighbor numbers and large series of bond angle powers
|
||||
void Fingerprint_bondscreened::do3bodyfeatureset_doubleneighborloop(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,ktype,kk,m,n;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int jtypes = atomtypes[1];
|
||||
int ktypes = atomtypes[2];
|
||||
int count=0;
|
||||
PairRANN::Simulation *sim = &pair->sims[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;
|
||||
// numneigh = sim->numneigh;
|
||||
// firstneigh = sim->firstneigh;
|
||||
i = ilist[ii];
|
||||
itype = pair->map[type[i]];
|
||||
int f = pair->net[itype].dimensions[0];
|
||||
// xtmp = x[i][0];
|
||||
// ytmp = x[i][1];
|
||||
// ztmp = x[i][2];
|
||||
// jlist = firstneigh[i];
|
||||
// jnum = numneigh[i];
|
||||
double expr[jnum][this->k];
|
||||
double y[jnum][3];
|
||||
double ri[jnum];
|
||||
double dfc[jnum];
|
||||
int kb = this->k;
|
||||
int mb = this->m;
|
||||
double Bijk[kb][mb];
|
||||
double c41[this->k];
|
||||
double c51[this->k];
|
||||
double c61[this->k];
|
||||
double ct[this->k];
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
if (Bij[jj]==false){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype){
|
||||
expr[jj][0]=0;
|
||||
continue;
|
||||
}
|
||||
// delx = xtmp - x[j][0];
|
||||
// dely = ytmp - x[j][1];
|
||||
// delz = ztmp - x[j][2];
|
||||
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("Neighbor list is invalid.");//usually results from nan somewhere.
|
||||
r1 = r1-trunc(r1);
|
||||
double *p0 = &expcuttable[(m1-1)*this->k];
|
||||
double *p1 = &expcuttable[m1*this->k];
|
||||
double *p2 = &expcuttable[(m1+1)*this->k];
|
||||
double *p3 = &expcuttable[(m1+2)*this->k];
|
||||
for (kk=0;kk<this->k;kk++){
|
||||
expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
|
||||
expr[jj][kk] *= Sik[jj];
|
||||
}
|
||||
double* q = &dfctable[m1-1];
|
||||
double* r2 = &rinvsqrttable[m1-1];
|
||||
dfc[jj] = 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])));
|
||||
ri[jj] = r2[1] + 0.5 * r1*(r2[2] - r2[0] + r1*(2.0*r2[0] - 5.0*r2[1] + 4.0*r2[2] - r2[3] + r1*(3.0*(r2[1] - r2[2]) + r2[3] - r2[0])));
|
||||
y[jj][0]=delx*ri[jj];
|
||||
y[jj][1]=dely*ri[jj];
|
||||
y[jj][2]=delz*ri[jj];
|
||||
}
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
if (Bij[jj]==false){continue;}
|
||||
if (expr[jj][0]==0)continue;
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype =tn[jj];
|
||||
if (jtypes != nelements && jtypes != jtype){
|
||||
continue;
|
||||
}
|
||||
for (n = 0;n<this->k;n++){
|
||||
ct[n] = alpha_k[n]/re;
|
||||
c41[n]=(-ct[n]+dfc[jj])*y[jj][0]+dSikx[jj];
|
||||
c51[n]=(-ct[n]+dfc[jj])*y[jj][1]+dSiky[jj];
|
||||
c61[n]=(-ct[n]+dfc[jj])*y[jj][2]+dSikz[jj];
|
||||
}
|
||||
for (n=0;n<kb;n++){for (m=0;m<mb;m++){Bijk[n][m]=0;}}
|
||||
// if (jtypes==ktypes){
|
||||
for (kk=0;kk< jnum; kk++){
|
||||
if (Bij[kk]==false){continue;}
|
||||
if (expr[kk][0]==0)continue;
|
||||
// int k1 = jlist[kk];
|
||||
// k1 &= NEIGHMASK;
|
||||
// ktype = pair->map[type[k1]];
|
||||
ktype = tn[kk];
|
||||
if (ktypes != nelements && ktypes != ktype){
|
||||
continue;
|
||||
}
|
||||
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<kb;n++){
|
||||
double dot1=expr[jj][n]*expr[kk][n];
|
||||
double c4 = c41[n];
|
||||
double c5 = c51[n];
|
||||
double c6 = c61[n];
|
||||
double ct2 = -ct[n]+dfc[kk];
|
||||
double c42 = ct2*y[kk][0]+dSikx[kk];
|
||||
double c52 = ct2*y[kk][1]+dSiky[kk];
|
||||
double c62 = ct2*y[kk][2]+dSikz[kk];
|
||||
//m=0
|
||||
Bijk[n][0]+=dot1;
|
||||
features[count]+=dot1;
|
||||
dfeaturesx[jj*f+count]+=dot1*c4;
|
||||
dfeaturesy[jj*f+count]+=dot1*c5;
|
||||
dfeaturesz[jj*f+count]+=dot1*c6;
|
||||
dfeaturesx[kk*f+count]+=dot1*c42;
|
||||
dfeaturesy[kk*f+count]+=dot1*c52;
|
||||
dfeaturesz[kk*f+count]+=dot1*c62;
|
||||
c4*=dot;
|
||||
c5*=dot;
|
||||
c6*=dot;
|
||||
c42*=dot;
|
||||
c52*=dot;
|
||||
c62*=dot;
|
||||
count++;
|
||||
for (m=1;m<mb;m++){
|
||||
double c7 = dot1*(m*c1+c4);
|
||||
double c8 = dot1*(m*c2+c5);
|
||||
double c9 = dot1*(m*c3+c6);
|
||||
dfeaturesx[jj*f+count]+=c7;
|
||||
dfeaturesy[jj*f+count]+=c8;
|
||||
dfeaturesz[jj*f+count]+=c9;
|
||||
dfeaturesx[kk*f+count]+=dot1*(m*c10+c42);
|
||||
dfeaturesy[kk*f+count]+=dot1*(m*c11+c52);
|
||||
dfeaturesz[kk*f+count]+=dot1*(m*c12+c62);
|
||||
dot1*=dot;
|
||||
features[count++]+=dot1;
|
||||
Bijk[n][m] += dot1;
|
||||
}
|
||||
}
|
||||
}
|
||||
// kk=jj;
|
||||
// if (ktypes == nelements || ktypes == jtype){
|
||||
// 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 c1 = 2*ri[jj]*y[kk][0]*(1-y[jj][0]*y[jj][0]);
|
||||
//// double c2 = 2*ri[jj]*y[kk][1]*(1-y[jj][1]*y[jj][1]);
|
||||
//// double c3 = 2*ri[jj]*y[kk][2]*(1-y[jj][2]*y[jj][2]);
|
||||
// for (n=0;n<kb;n++){
|
||||
// double dot1=expr[jj][n]*expr[kk][n]*pair->Sik[kk]*pair->Sik[jj];
|
||||
// double c4 = c41[n]+4*pair->dSikx[kk];
|
||||
// double c5 = c51[n]+4*pair->dSiky[kk];
|
||||
// double c6 = c61[n]+4*pair->dSikz[kk];
|
||||
// //m=0
|
||||
// features[count]+=dot1;
|
||||
// Bijk[n][0]+=dot1;
|
||||
// dfeaturesx[jj*f+count]+=dot1*c4;
|
||||
// dfeaturesy[jj*f+count]+=dot1*c5;
|
||||
// dfeaturesz[jj*f+count]+=dot1*c6;
|
||||
// c4*=dot;
|
||||
// c5*=dot;
|
||||
// c6*=dot;
|
||||
// count++;
|
||||
// for (m=1;m<mb;m++){
|
||||
// double c7 = dot1*(m*c1+c4);
|
||||
// double c8 = dot1*(m*c2+c5);
|
||||
// double c9 = dot1*(m*c3+c6);
|
||||
// dfeaturesx[jj*f+count]+=c7;
|
||||
// dfeaturesy[jj*f+count]+=c8;
|
||||
// dfeaturesz[jj*f+count]+=c9;
|
||||
// dot1*=dot;
|
||||
// features[count++]+=dot1;
|
||||
// Bijk[n][m] += dot1;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
for (kk=0;kk<jnum;kk++){
|
||||
if (Bij[kk]==false){continue;}
|
||||
if (expr[kk][0]==0)continue;
|
||||
// int k1 = jlist[kk];
|
||||
// k1 &= NEIGHMASK;
|
||||
// ktype = pair->map[type[k1]];
|
||||
ktype = tn[kk];
|
||||
if (ktypes != nelements && ktypes != ktype){
|
||||
continue;
|
||||
}
|
||||
count = startingneuron;
|
||||
for (n=0;n<kb;n++){
|
||||
for (m=0;m<mb;m++){
|
||||
dfeaturesx[kk*f+count]+=2*Bijk[n][m]*dSijkx[jj*jnum+kk];
|
||||
dfeaturesy[kk*f+count]+=2*Bijk[n][m]*dSijky[jj*jnum+kk];
|
||||
dfeaturesz[kk*f+count]+=2*Bijk[n][m]*dSijkz[jj*jnum+kk];
|
||||
count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
// }
|
||||
// else {
|
||||
// for (kk=0;kk<jnum; kk++){
|
||||
// if (pair->Bij[kk]==false){continue;}
|
||||
// if (expr[kk][0]==0)continue;
|
||||
// int k1 = jlist[kk];
|
||||
// k1 &= NEIGHMASK;
|
||||
// ktype = pair->map[type[k1]];
|
||||
// if (ktypes != nelements && ktypes != ktype){
|
||||
// continue;
|
||||
// }
|
||||
// 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<kb;n++){
|
||||
// double dot1=expr[jj][n]*expr[kk][n]*pair->Sik[kk]*pair->Sik[jj];
|
||||
// double c4 = c41[n]/2;
|
||||
// double c5 = c51[n]/2;
|
||||
// double c6 = c61[n]/2;
|
||||
// double ct2 = -ct[n]/2+dfc[kk];
|
||||
// double c42 = ct2*y[kk][0]+pair->dSikx[kk];
|
||||
// double c52 = ct2*y[kk][1]+pair->dSiky[kk];
|
||||
// double c62 = ct2*y[kk][2]+pair->dSikz[kk];
|
||||
// //m=0
|
||||
// features[count]+=dot1;
|
||||
// dfeaturesx[jj*f+count]+=dot1*c4;
|
||||
// dfeaturesy[jj*f+count]+=dot1*c5;
|
||||
// dfeaturesz[jj*f+count]+=dot1*c6;
|
||||
// dfeaturesx[kk*f+count]+=dot1*c42;
|
||||
// dfeaturesy[kk*f+count]+=dot1*c52;
|
||||
// dfeaturesz[kk*f+count]+=dot1*c62;
|
||||
// c4*=dot;
|
||||
// c5*=dot;
|
||||
// c6*=dot;
|
||||
// c42*=dot;
|
||||
// c52*=dot;
|
||||
// c62*=dot;
|
||||
// count++;
|
||||
// for (m=1;m<mb;m++){
|
||||
// double c7 = dot1*(m*c1+c4);
|
||||
// double c8 = dot1*(m*c2+c5);
|
||||
// double c9 = dot1*(m*c3+c6);
|
||||
// dfeaturesx[jj*f+count]+=c7;
|
||||
// dfeaturesy[jj*f+count]+=c8;
|
||||
// dfeaturesz[jj*f+count]+=c9;
|
||||
// dfeaturesx[kk*f+count]+=dot1*(m*c10+c42);
|
||||
// dfeaturesy[kk*f+count]+=dot1*(m*c11+c52);
|
||||
// dfeaturesz[kk*f+count]+=dot1*(m*c12+c62);
|
||||
// dot1*=dot;
|
||||
// features[count++]+=dot1;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
}
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (Bij[jj]==false){continue;}
|
||||
if (expr[jj][0]==0){continue;}
|
||||
count = startingneuron;
|
||||
for (n=0;n<kb;n++){
|
||||
for (m=0;m<mb;m++){
|
||||
dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
|
||||
dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
|
||||
dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
|
||||
count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int Fingerprint_bondscreened::factorial(int n) {
|
||||
if ((n==0)||(n==1))
|
||||
return 1;
|
||||
else
|
||||
return n*factorial(n-1);
|
||||
}
|
||||
67
src/fingerprint_bondscreened.h
Normal file
67
src/fingerprint_bondscreened.h
Normal file
@ -0,0 +1,67 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#ifdef FINGERPRINT_CLASS
|
||||
|
||||
FingerprintStyle(bondscreened,Fingerprint_bondscreened)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef FINGERPRINT_BONDSCREENED_H_
|
||||
#define FINGERPRINT_BONDSCREENED_H_
|
||||
|
||||
#include "fingerprint.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Fingerprint_bondscreened : public Fingerprint {
|
||||
public:
|
||||
Fingerprint_bondscreened(PairRANN *);
|
||||
~Fingerprint_bondscreened();
|
||||
bool parse_values(char*,char*);
|
||||
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*);
|
||||
int factorial(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 k;
|
||||
int m;
|
||||
int **Mf;
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
#endif /* FINGERPRINT_BOND_H_ */
|
||||
1019
src/fingerprint_bondscreenedspin.cpp
Normal file
1019
src/fingerprint_bondscreenedspin.cpp
Normal file
File diff suppressed because it is too large
Load Diff
67
src/fingerprint_bondscreenedspin.h
Normal file
67
src/fingerprint_bondscreenedspin.h
Normal file
@ -0,0 +1,67 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#ifdef FINGERPRINT_CLASS
|
||||
|
||||
FingerprintStyle(bondscreenedspin,Fingerprint_bondscreenedspin)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef FINGERPRINT_BONDSCREENEDSPIN_H_
|
||||
#define FINGERPRINT_BONDSCREENEDSPIN_H_
|
||||
|
||||
#include "fingerprint.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Fingerprint_bondscreenedspin : public Fingerprint {
|
||||
public:
|
||||
Fingerprint_bondscreenedspin(PairRANN *);
|
||||
~Fingerprint_bondscreenedspin();
|
||||
bool parse_values(char*,char*);
|
||||
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*);
|
||||
int factorial(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 k;
|
||||
int m;
|
||||
int **Mf;
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
#endif /* FINGERPRINT_BOND_H_ */
|
||||
1015
src/fingerprint_bondspin.cpp
Normal file
1015
src/fingerprint_bondspin.cpp
Normal file
File diff suppressed because it is too large
Load Diff
67
src/fingerprint_bondspin.h
Normal file
67
src/fingerprint_bondspin.h
Normal file
@ -0,0 +1,67 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#ifdef FINGERPRINT_CLASS
|
||||
|
||||
FingerprintStyle(bondspin,Fingerprint_bondspin)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef FINGERPRINT_BONDSPIN_H_
|
||||
#define FINGERPRINT_BONDSPIN_H_
|
||||
|
||||
#include "fingerprint.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Fingerprint_bondspin : public Fingerprint {
|
||||
public:
|
||||
Fingerprint_bondspin(PairRANN *);
|
||||
~Fingerprint_bondspin();
|
||||
bool parse_values(char*,char*);
|
||||
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*);
|
||||
int factorial(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 k;
|
||||
int m;
|
||||
int **Mf;
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
#endif /* FINGERPRINT_BOND_H_ */
|
||||
253
src/fingerprint_radial.cpp
Normal file
253
src/fingerprint_radial.cpp
Normal file
@ -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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
|
||||
#include "fingerprint_radial.h"
|
||||
#include "fingerprint.h"
|
||||
#include <math.h>
|
||||
#include <cmath>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <iostream>
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
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;
|
||||
n = 0;
|
||||
o = 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(char* constant,char * line1){
|
||||
int l;
|
||||
char **words=new char *[MAXLINE];
|
||||
int nwords;
|
||||
nwords=0;
|
||||
words[nwords++] = strtok(line1,": ,\t\n");
|
||||
while ((words[nwords++] = strtok(NULL,": ,\t\n"))) continue;
|
||||
nwords -= 1;
|
||||
if (strcmp(constant,"re")==0){
|
||||
re = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"rc")==0){
|
||||
rc = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"alpha")==0){
|
||||
delete [] alpha;
|
||||
alpha = new double [nwords];
|
||||
for (l=0;l<nwords;l++){
|
||||
alpha[l]=strtod(words[l],NULL);
|
||||
}
|
||||
}
|
||||
else if (strcmp(constant,"dr")==0){
|
||||
dr = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"n")==0){
|
||||
n = strtol(line1,NULL,10);
|
||||
}
|
||||
else if (strcmp(constant,"o")==0){
|
||||
o = strtol(line1,NULL,10);
|
||||
}
|
||||
else pair->errorf("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.
|
||||
delete [] words;
|
||||
if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && n!=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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:alpha:\n",style,id);
|
||||
for (i=0;i<(n-o+1);i++){
|
||||
fprintf(fid,"%f ",alpha[i]);
|
||||
}
|
||||
fprintf(fid,"\n");
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:o:\n",style,id);
|
||||
fprintf(fid,"%d\n",o);
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:n:\n",style,id);
|
||||
fprintf(fid,"%d\n",n);
|
||||
}
|
||||
|
||||
//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<=(n-o);m++){
|
||||
radialtable[k*(n-o+1)+m]=pow(sqrt(r1)/re,m+o)*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;j<n_body_type;j++){atomtypes[j] = i[j];}
|
||||
this->id = id;
|
||||
}
|
||||
|
||||
void Fingerprint_radial::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 nelements = pair->nelements;
|
||||
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;
|
||||
// numneigh = sim->numneigh;
|
||||
// firstneigh = sim->firstneigh;
|
||||
// xtmp = x[i][0];
|
||||
// ytmp = x[i][1];
|
||||
// ztmp = x[i][2];
|
||||
// jlist = firstneigh[i];
|
||||
// jnum = numneigh[i];
|
||||
//loop over neighbors
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype =tn[jj];
|
||||
if (this->atomtypes[1] != nelements && this->atomtypes[1] != jtype)continue;
|
||||
// delx = xtmp - x[j][0];
|
||||
// dely = ytmp - x[j][1];
|
||||
// delz = ztmp - x[j][2];
|
||||
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("invalid neighbor radius!");}
|
||||
if (radialtable[m1]==0){continue;}
|
||||
//cubic interpolation from tables
|
||||
double *p1 = &radialtable[m1*(n-o+1)];
|
||||
double *p2 = &radialtable[(m1+1)*(n-o+1)];
|
||||
double *p3 = &radialtable[(m1+2)*(n-o+1)];
|
||||
double *p0 = &radialtable[(m1-1)*(n-o+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<=(n-o);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+o)/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 n-o+1;
|
||||
}
|
||||
56
src/fingerprint_radial.h
Normal file
56
src/fingerprint_radial.h
Normal file
@ -0,0 +1,56 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#ifdef FINGERPRINT_CLASS
|
||||
|
||||
FingerprintStyle(radial,Fingerprint_radial)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef FINGERPRINT_RADIAL_H_
|
||||
#define FINGERPRINT_RADIAL_H_
|
||||
|
||||
#include "fingerprint.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Fingerprint_radial : public Fingerprint {
|
||||
public:
|
||||
Fingerprint_radial(PairRANN *);
|
||||
~Fingerprint_radial();
|
||||
bool parse_values(char*,char*);
|
||||
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 n;//highest term
|
||||
int o;//lowest term
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif /* FINGERPRINT_RADIAL_H_ */
|
||||
266
src/fingerprint_radialscreened.cpp
Normal file
266
src/fingerprint_radialscreened.cpp
Normal file
@ -0,0 +1,266 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
|
||||
#include "fingerprint_radialscreened.h"
|
||||
#include <cmath>
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <iostream>
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
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;
|
||||
n = 0;
|
||||
o = 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(char* constant,char * line1){
|
||||
int l;
|
||||
char **words=new char *[MAXLINE];
|
||||
int nwords;
|
||||
nwords=0;
|
||||
words[nwords++] = strtok(line1,": ,\t\n");
|
||||
while ((words[nwords++] = strtok(NULL,": ,\t\n"))) continue;
|
||||
nwords -= 1;
|
||||
if (strcmp(constant,"re")==0){
|
||||
re = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"rc")==0){
|
||||
rc = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"alpha")==0){
|
||||
delete [] alpha;
|
||||
alpha = new double [nwords];
|
||||
for (l=0;l<nwords;l++){
|
||||
alpha[l]=strtod(words[l],NULL);
|
||||
}
|
||||
}
|
||||
else if (strcmp(constant,"dr")==0){
|
||||
dr = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"n")==0){
|
||||
n = strtol(line1,NULL,10);
|
||||
}
|
||||
else if (strcmp(constant,"o")==0){
|
||||
o = strtol(line1,NULL,10);
|
||||
}
|
||||
else pair->errorf("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.
|
||||
delete [] words;
|
||||
if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && n!=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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:alpha:\n",style,id);
|
||||
for (i=0;i<(n-o+1);i++){
|
||||
fprintf(fid,"%f ",alpha[i]);
|
||||
}
|
||||
fprintf(fid,"\n");
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:o:\n",style,id);
|
||||
fprintf(fid,"%d\n",o);
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:n:\n",style,id);
|
||||
fprintf(fid,"%d\n",n);
|
||||
}
|
||||
|
||||
//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<=(n-o);m++){
|
||||
radialtable[k*(n-o+1)+m]=pow(sqrt(r1)/re,m+o)*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;j<n_body_type;j++){atomtypes[j] = i[j];}
|
||||
this->id = id;
|
||||
}
|
||||
|
||||
void Fingerprint_radialscreened::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 nelements = pair->nelements;
|
||||
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;
|
||||
// numneigh = sim->numneigh;
|
||||
// firstneigh = sim->firstneigh;
|
||||
// xtmp = x[i][0];
|
||||
// ytmp = x[i][1];
|
||||
// ztmp = x[i][2];
|
||||
// jlist = firstneigh[i];
|
||||
// jnum = numneigh[i];
|
||||
//loop over neighbors
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
if (Bij[jj]==false){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (this->atomtypes[1] != nelements && this->atomtypes[1] != jtype)continue;
|
||||
// delx = xtmp - x[j][0];
|
||||
// dely = ytmp - x[j][1];
|
||||
// delz = ztmp - x[j][2];
|
||||
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("invalid neighbor radius!");}
|
||||
if (radialtable[m1]==0){continue;}
|
||||
//cubic interpolation from tables
|
||||
double *p1 = &radialtable[m1*(n-o+1)];
|
||||
double *p2 = &radialtable[(m1+1)*(n-o+1)];
|
||||
double *p3 = &radialtable[(m1+2)*(n-o+1)];
|
||||
double *p0 = &radialtable[(m1-1)*(n-o+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<=(n-o);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+o)/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<jnum;kk++){
|
||||
if (Bij[kk]==false){continue;}
|
||||
dfeaturesx[kk*f+count]+=rt*dSijkx[jj*jnum+kk];
|
||||
dfeaturesy[kk*f+count]+=rt*dSijky[jj*jnum+kk];
|
||||
dfeaturesz[kk*f+count]+=rt*dSijkz[jj*jnum+kk];
|
||||
}
|
||||
count++;
|
||||
}
|
||||
}
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (Bij[jj]==false){continue;}
|
||||
count = startingneuron;
|
||||
for (l=0;l<=(n-o);l++){
|
||||
dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
|
||||
dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
|
||||
dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
int Fingerprint_radialscreened::get_length()
|
||||
{
|
||||
return n-o+1;
|
||||
}
|
||||
59
src/fingerprint_radialscreened.h
Normal file
59
src/fingerprint_radialscreened.h
Normal file
@ -0,0 +1,59 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#ifdef FINGERPRINT_CLASS
|
||||
|
||||
FingerprintStyle(radialscreened,Fingerprint_radialscreened)
|
||||
|
||||
#else
|
||||
#ifndef FINGERPRINT_RADIALSCREENED_H_
|
||||
#define FINGERPRINT_RADIALSCREENED_H_
|
||||
|
||||
|
||||
#include "fingerprint.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Fingerprint_radialscreened : public Fingerprint {
|
||||
public:
|
||||
Fingerprint_radialscreened(PairRANN *);
|
||||
~Fingerprint_radialscreened();
|
||||
bool parse_values(char*,char*);
|
||||
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 n;//highest term
|
||||
int o;//lowest term
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
#endif /* FINGERPRINT_RADIALSCREENED_H_ */
|
||||
279
src/fingerprint_radialscreenedspin.cpp
Normal file
279
src/fingerprint_radialscreenedspin.cpp
Normal file
@ -0,0 +1,279 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
|
||||
#include "fingerprint_radialscreenedspin.h"
|
||||
#include <cmath>
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <iostream>
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
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;
|
||||
n = 0;
|
||||
o = 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(char* constant,char * line1){
|
||||
int l;
|
||||
char **words=new char *[MAXLINE];
|
||||
int nwords;
|
||||
nwords=0;
|
||||
words[nwords++] = strtok(line1,": ,\t\n");
|
||||
while ((words[nwords++] = strtok(NULL,": ,\t\n"))) continue;
|
||||
nwords -= 1;
|
||||
if (strcmp(constant,"re")==0){
|
||||
re = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"rc")==0){
|
||||
rc = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"alpha")==0){
|
||||
delete [] alpha;
|
||||
alpha = new double [nwords];
|
||||
for (l=0;l<nwords;l++){
|
||||
alpha[l]=strtod(words[l],NULL);
|
||||
}
|
||||
}
|
||||
else if (strcmp(constant,"dr")==0){
|
||||
dr = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"n")==0){
|
||||
n = strtol(line1,NULL,10);
|
||||
}
|
||||
else if (strcmp(constant,"o")==0){
|
||||
o = strtol(line1,NULL,10);
|
||||
}
|
||||
else pair->errorf("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.
|
||||
delete [] words;
|
||||
if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && n!=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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:alpha:\n",style,id);
|
||||
for (i=0;i<(n-o+1);i++){
|
||||
fprintf(fid,"%f ",alpha[i]);
|
||||
}
|
||||
fprintf(fid,"\n");
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:o:\n",style,id);
|
||||
fprintf(fid,"%d\n",o);
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:n:\n",style,id);
|
||||
fprintf(fid,"%d\n",n);
|
||||
}
|
||||
|
||||
//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<=(n-o);m++){
|
||||
radialtable[k*(n-o+1)+m]=pow(sqrt(r1)/re,m+o)*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;j<n_body_type;j++){atomtypes[j] = i[j];}
|
||||
this->id = id;
|
||||
}
|
||||
|
||||
void Fingerprint_radialscreenedspin::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 nelements = pair->nelements;
|
||||
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];
|
||||
// numneigh = sim->numneigh;
|
||||
// firstneigh = sim->firstneigh;
|
||||
// xtmp = x[i][0];
|
||||
// ytmp = x[i][1];
|
||||
// ztmp = x[i][2];
|
||||
// jlist = firstneigh[i];
|
||||
// jnum = numneigh[i];
|
||||
//loop over neighbors
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
if (Bij[jj]==false){continue;}
|
||||
// j = jlist[jj];
|
||||
// j &= NEIGHMASK;
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype = tn[jj];
|
||||
if (this->atomtypes[1] != nelements && this->atomtypes[1] != jtype)continue;
|
||||
// delx = xtmp - x[j][0];
|
||||
// dely = ytmp - x[j][1];
|
||||
// delz = ztmp - x[j][2];
|
||||
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("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*(n-o+1)];
|
||||
double *p2 = &radialtable[(m1+1)*(n-o+1)];
|
||||
double *p3 = &radialtable[(m1+2)*(n-o+1)];
|
||||
double *p0 = &radialtable[(m1-1)*(n-o+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<=(n-o);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+o)/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<jnum;kk++){
|
||||
if (Bij[kk]==false){continue;}
|
||||
dfeaturesx[kk*f+count]+=rt*dSijkx[jj*jnum+kk];
|
||||
dfeaturesy[kk*f+count]+=rt*dSijky[jj*jnum+kk];
|
||||
dfeaturesz[kk*f+count]+=rt*dSijkz[jj*jnum+kk];
|
||||
}
|
||||
count++;
|
||||
}
|
||||
}
|
||||
for (jj=0;jj<jnum;jj++){
|
||||
if (Bij[jj]==false){continue;}
|
||||
count = startingneuron;
|
||||
for (l=0;l<=(n-o);l++){
|
||||
dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
|
||||
dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
|
||||
dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
int Fingerprint_radialscreenedspin::get_length()
|
||||
{
|
||||
return n-o+1;
|
||||
}
|
||||
59
src/fingerprint_radialscreenedspin.h
Normal file
59
src/fingerprint_radialscreenedspin.h
Normal file
@ -0,0 +1,59 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#ifdef FINGERPRINT_CLASS
|
||||
|
||||
FingerprintStyle(radialscreenedspin,Fingerprint_radialscreenedspin)
|
||||
|
||||
#else
|
||||
#ifndef FINGERPRINT_RADIALSCREENEDSPIN_H_
|
||||
#define FINGERPRINT_RADIALSCREENEDSPIN_H_
|
||||
|
||||
|
||||
#include "fingerprint.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Fingerprint_radialscreenedspin : public Fingerprint {
|
||||
public:
|
||||
Fingerprint_radialscreenedspin(PairRANN *);
|
||||
~Fingerprint_radialscreenedspin();
|
||||
bool parse_values(char*,char*);
|
||||
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 n;//highest term
|
||||
int o;//lowest term
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
#endif /* FINGERPRINT_RADIALSCREENED_H_ */
|
||||
265
src/fingerprint_radialspin.cpp
Normal file
265
src/fingerprint_radialspin.cpp
Normal file
@ -0,0 +1,265 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
|
||||
#include "fingerprint_radialspin.h"
|
||||
|
||||
#include "fingerprint.h"
|
||||
#include <math.h>
|
||||
#include <cmath>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <iostream>
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
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;
|
||||
n = 0;
|
||||
o = 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(char* constant,char * line1){
|
||||
int l;
|
||||
char **words=new char *[MAXLINE];
|
||||
int nwords;
|
||||
nwords=0;
|
||||
words[nwords++] = strtok(line1,": ,\t\n");
|
||||
while ((words[nwords++] = strtok(NULL,": ,\t\n"))) continue;
|
||||
nwords -= 1;
|
||||
if (strcmp(constant,"re")==0){
|
||||
re = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"rc")==0){
|
||||
rc = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"alpha")==0){
|
||||
delete [] alpha;
|
||||
alpha = new double [nwords];
|
||||
for (l=0;l<nwords;l++){
|
||||
alpha[l]=strtod(words[l],NULL);
|
||||
}
|
||||
}
|
||||
else if (strcmp(constant,"dr")==0){
|
||||
dr = strtod(line1,NULL);
|
||||
}
|
||||
else if (strcmp(constant,"n")==0){
|
||||
n = strtol(line1,NULL,10);
|
||||
}
|
||||
else if (strcmp(constant,"o")==0){
|
||||
o = strtol(line1,NULL,10);
|
||||
}
|
||||
else pair->errorf("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.
|
||||
delete [] words;
|
||||
if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && n!=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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:alpha:\n",style,id);
|
||||
for (i=0;i<(n-o+1);i++){
|
||||
fprintf(fid,"%f ",alpha[i]);
|
||||
}
|
||||
fprintf(fid,"\n");
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[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;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:o:\n",style,id);
|
||||
fprintf(fid,"%d\n",o);
|
||||
fprintf(fid,"fingerprintconstants:");
|
||||
fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
|
||||
for (i=1;i<n_body_type;i++){
|
||||
fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
|
||||
}
|
||||
fprintf(fid,":%s_%d:n:\n",style,id);
|
||||
fprintf(fid,"%d\n",n);
|
||||
}
|
||||
|
||||
//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<=(n-o);m++){
|
||||
radialtable[k*(n-o+1)+m]=pow(sqrt(r1)/re,m+o)*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;j<n_body_type;j++){atomtypes[j] = i[j];}
|
||||
this->id = id;
|
||||
}
|
||||
|
||||
void Fingerprint_radialspin::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 nelements = pair->nelements;
|
||||
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;
|
||||
double *si = sim->s[i];
|
||||
// numneigh = sim->numneigh;
|
||||
firstneigh = sim->firstneigh;
|
||||
// xtmp = x[i][0];
|
||||
// ytmp = x[i][1];
|
||||
// ztmp = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
// jnum = numneigh[i];
|
||||
//loop over neighbors
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jl[jj];
|
||||
// jtype = pair->map[type[j]];
|
||||
jtype =tn[jj];
|
||||
if (this->atomtypes[1] != nelements && this->atomtypes[1] != jtype)continue;
|
||||
// delx = xtmp - x[j][0];
|
||||
// dely = ytmp - x[j][1];
|
||||
// delz = ztmp - x[j][2];
|
||||
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("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*(n-o+1)];
|
||||
double *p2 = &radialtable[(m1+1)*(n-o+1)];
|
||||
double *p3 = &radialtable[(m1+2)*(n-o+1)];
|
||||
double *p0 = &radialtable[(m1-1)*(n-o+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<=(n-o);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+o)/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 n-o+1;
|
||||
}
|
||||
56
src/fingerprint_radialspin.h
Normal file
56
src/fingerprint_radialspin.h
Normal file
@ -0,0 +1,56 @@
|
||||
/* -*- 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 author: Christopher Barrett (MSU) barrett@me.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
#ifdef FINGERPRINT_CLASS
|
||||
|
||||
FingerprintStyle(radialspin,Fingerprint_radialspin)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef FINGERPRINT_RADIALspin_H_
|
||||
#define FINGERPRINT_RADIALspin_H_
|
||||
|
||||
#include "fingerprint.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class Fingerprint_radialspin : public Fingerprint {
|
||||
public:
|
||||
Fingerprint_radialspin(PairRANN *);
|
||||
~Fingerprint_radialspin();
|
||||
bool parse_values(char*,char*);
|
||||
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 n;//highest term
|
||||
int o;//lowest term
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif /* FINGERPRINT_RADIAL_H_ */
|
||||
2423
src/pair_rann.cpp
Normal file
2423
src/pair_rann.cpp
Normal file
File diff suppressed because it is too large
Load Diff
153
src/pair_rann.h
Normal file
153
src/pair_rann.h
Normal file
@ -0,0 +1,153 @@
|
||||
/* -*- 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@cavs.msstate.edu
|
||||
----------------------------------------------------------------------*/
|
||||
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(rann,PairRANN)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_RANN
|
||||
#define LMP_PAIR_RANN
|
||||
|
||||
#define MAXLINE 4096
|
||||
|
||||
#include "neigh_list.h"
|
||||
#include "pair.h"
|
||||
#include <map>
|
||||
#include <math.h>
|
||||
#include <string>
|
||||
#include <sys/resource.h>
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
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(char*);
|
||||
int count_words(char *);
|
||||
//black magic for modular fingerprints and activations
|
||||
class Activation ***activation;
|
||||
class Fingerprint ***fingerprints;
|
||||
typedef Fingerprint *(*FingerprintCreator)(PairRANN *);
|
||||
typedef Activation *(*ActivationCreator)(PairRANN *);
|
||||
typedef std::map<std::string,FingerprintCreator> FingerprintCreatorMap;
|
||||
typedef std::map<std::string,ActivationCreator> ActivationCreatorMap;
|
||||
FingerprintCreatorMap *fingerprint_map;
|
||||
ActivationCreatorMap *activation_map;
|
||||
Fingerprint * create_fingerprint(const char *);
|
||||
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 *Sik,*dSikx,*dSiky,*dSikz,*dSijkx,*dSijky,*dSijkz;
|
||||
// bool *Bij;
|
||||
double *screening_min;
|
||||
double *screening_max;
|
||||
bool **weightdefined;
|
||||
bool **biasdefined;
|
||||
|
||||
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;
|
||||
};
|
||||
Simulation *sims;
|
||||
|
||||
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)
|
||||
};
|
||||
NNarchitecture *net;//array of networks, 1 for each element.
|
||||
|
||||
private:
|
||||
template <typename T> static Fingerprint *fingerprint_creator(PairRANN *);
|
||||
template <typename T> static Activation *activation_creator(PairRANN *);
|
||||
//new functions
|
||||
void allocate(char **);//called after reading element list, but before reading the rest of the potential
|
||||
void read_file(char *);//read potential file
|
||||
void read_atom_types(char **,char *);
|
||||
void read_mass(char **,char *);
|
||||
void read_fpe(char**,char *);//fingerprints per element. Count total fingerprints defined for each 1st element in element combinations
|
||||
void read_fingerprints(char **,int,char *);
|
||||
void read_fingerprint_constants(char **,int,char *);
|
||||
void read_network_layers(char**,char*);//include input and output layer (hidden layers + 2)
|
||||
void read_layer_size(char**,char*);
|
||||
void read_weight(char**,char*,FILE*);//weights should be formatted as properly shaped matrices
|
||||
void read_bias(char**,char*,FILE*);//biases should be formatted as properly shaped vectors
|
||||
void read_activation_functions(char**,char*);
|
||||
void read_screening(char**,int, char*);
|
||||
bool check_potential();//after finishing reading potential file
|
||||
void propagateforward(double *,double *,double *,double *,double *,double **,double **,int,int,int*);//called by compute to get force and energy
|
||||
void propagateforwardspin(double *,double *,double *,double *,double *,double *,double *,double *,double **,double **,double**,int,int,int*);//called by compute to get force and energy
|
||||
void screen(double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int *,int);
|
||||
void testdfeatures(double *,double *, double *, double *,double *,double *, double *,int);
|
||||
void testdenergy();
|
||||
void cull_neighbor_list(double *,double *,double *,int *,int *,int *,int,int);
|
||||
void screen_neighbor_list(double *,double *,double *,int *,int *,int *,int,int,bool*,double*,double*,double*,double*,double*,double*,double*);
|
||||
void update_stack_size();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
2
src/style_activation.h
Normal file
2
src/style_activation.h
Normal file
@ -0,0 +1,2 @@
|
||||
#include "activation_linear.h"
|
||||
#include "activation_sigI.h"
|
||||
8
src/style_fingerprint.h
Normal file
8
src/style_fingerprint.h
Normal file
@ -0,0 +1,8 @@
|
||||
#include "fingerprint_bond.h"
|
||||
#include "fingerprint_bondscreened.h"
|
||||
#include "fingerprint_bondscreenedspin.h"
|
||||
#include "fingerprint_bondspin.h"
|
||||
#include "fingerprint_radial.h"
|
||||
#include "fingerprint_radialscreened.h"
|
||||
#include "fingerprint_radialscreenedspin.h"
|
||||
#include "fingerprint_radialspin.h"
|
||||
Reference in New Issue
Block a user