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

refactor to use TabularFunction class shared with pair style bop

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-05-04 22:26:56 -04:00
parent 8285e71fd9
commit cdc3e2bad9
2 changed files with 196 additions and 397 deletions
Download Patch File Download Diff File Expand all files Collapse all files

330
src/MANYBODY/pair_polymorphic.cpp
View File

@ -16,6 +16,9 @@
This modifies from pair_tersoff.cpp by Aidan Thompson (SNL)
------------------------------------------------------------------------- */
// uncomment define to enable writing table files for debugging
// #define LMP_POLYMORPHIC_WRITE_TABLES 1
#include "pair_polymorphic.h"
#include "atom.h"
@ -28,6 +31,7 @@
#include "neigh_request.h"
#include "neighbor.h"
#include "potential_file_reader.h"
#include "tabular_function.h"
#include "tokenizer.h"
#include <cmath>
@ -39,6 +43,42 @@ using namespace MathExtra;
#define MAXLINE 1024
#define DELTA 4
/* ---------------------------------------------------------------------- */
PairPolymorphic::PairParameters::PairParameters()
{
cut = 0.0;
cutsq = 0.0;
xi = 0.0;
U = nullptr;
V = nullptr;
W = nullptr;
F = nullptr;
}
PairPolymorphic::PairParameters::~PairParameters()
{
delete U;
delete V;
delete W;
delete F;
}
/* ---------------------------------------------------------------------- */
PairPolymorphic::TripletParameters::TripletParameters()
{
P = nullptr;
G = nullptr;
}
PairPolymorphic::TripletParameters::~TripletParameters()
{
delete P;
delete G;
}
/* ---------------------------------------------------------------------- */
PairPolymorphic::PairPolymorphic(LAMMPS *lmp) : Pair(lmp)
@ -242,7 +282,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
}
if (rsq >= (p.U)->get_xmaxsq() || (p.U)->get_vmax() <= epsilon) continue;
(p.U)->value(r0,evdwl,eflag,fpair,1);
(p.U)->value2(r0,evdwl,eflag,fpair,1);
fpair = -fpair/r0;
f[i][0] += delx*fpair;
@ -274,14 +314,14 @@ void PairPolymorphic::compute(int eflag, int vflag)
iparam_kk = elem2param[ktype][ktype];
PairParameters & q = pairParameters[iparam_kk];
(q.W)->value(drW[kk],wfac,1,fpair,0);
(q.W)->value2(drW[kk],wfac,1,fpair,0);
zeta_ij += wfac;
}
// pairwise force due to zeta
(p.F)->value(zeta_ij,bij,1,bij_d,1);
(p.F)->value2(zeta_ij,bij,1,bij_d,1);
prefactor = 0.5* bij_d;
if (eflag) evdwl = -0.5*bij;
@ -301,7 +341,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
iparam_kk = elem2param[ktype][ktype];
PairParameters & q = pairParameters[iparam_kk];
(q.W)->value(drW[kk],wfac,0,fpair,1);
(q.W)->value2(drW[kk],wfac,0,fpair,1);
fpair = -prefactor*fpair/drW[kk];
f[i][0] += delr2[0]*fpair;
@ -356,17 +396,17 @@ void PairPolymorphic::compute(int eflag, int vflag)
iparam_ik = elem2param[itype][ktype];
PairParameters & q = pairParameters[iparam_ik];
(q.W)->value(r2,wfac,1,fpair,0);
(trip.P)->value(r1-(p.xi)*r2,pfac,1,fpair,0);
(trip.G)->value(costheta,gfac,1,fpair,0);
(q.W)->value2(r2,wfac,1,fpair,0);
(trip.P)->value2(r1-(p.xi)*r2,pfac,1,fpair,0);
(trip.G)->value2(costheta,gfac,1,fpair,0);
zeta_ij += wfac*pfac*gfac;
}
// pairwise force due to zeta
(p.V)->value(r1,fa,1,fa_d,1);
(p.F)->value(zeta_ij,bij,1,bij_d,1);
(p.V)->value2(r1,fa,1,fa_d,1);
(p.F)->value2(zeta_ij,bij,1,bij_d,1);
fpair = -0.5*bij*fa_d / r1;
prefactor = 0.5* fa * bij_d;
if (eflag) evdwl = -0.5*bij*fa;
@ -606,8 +646,8 @@ void PairPolymorphic::read_file(char *file)
reader->next_dvector(singletable, nr);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.U = new tabularFunction(nr,0.0,p.cut);
(p.U)->set_values(nr,0.0,p.cut,singletable,epsilon);
p.U = new TabularFunction;
(p.U)->set_values(nr,0.0,p.cut,singletable);
}
for (int i = 0; i < npair; i++) { // V
PairParameters & p = pairParameters[i];
@ -615,8 +655,8 @@ void PairPolymorphic::read_file(char *file)
reader->next_dvector(singletable, nr);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.V = new tabularFunction(nr,0.0,p.cut);
(p.V)->set_values(nr,0.0,p.cut,singletable,epsilon);
p.V = new TabularFunction;
(p.V)->set_values(nr,0.0,p.cut,singletable);
}
for (int i = 0; i < npair; i++) { // W
PairParameters & p = pairParameters[i];
@ -624,8 +664,8 @@ void PairPolymorphic::read_file(char *file)
reader->next_dvector(singletable, nr);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.W = new tabularFunction(nr,0.0,p.cut);
(p.W)->set_values(nr,0.0,p.cut,singletable,epsilon);
p.W = new TabularFunction;
(p.W)->set_values(nr,0.0,p.cut,singletable);
}
cutmax = 0.0;
@ -643,25 +683,25 @@ void PairPolymorphic::read_file(char *file)
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
for (int i = 0; i < nelements; i++) {
TripletParameters & p = tripletParameters[i*nelements*nelements+j*nelements+j];
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
p.P = new TabularFunction;
(p.P)->set_values(nr,-cutmax,cutmax,singletable);
}
}
for (int j = 0; j < nelements-1; j++) { // P
for (int k = j+1; k < nelements; k++) {
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
for (int k = j+1; k < nelements; k++) {
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
for (int i = 0; i < nelements; i++) {
TripletParameters & p = tripletParameters[i*nelements*nelements+j*nelements+k];
p.P = new TabularFunction;
(p.P)->set_values(nr,-cutmax,cutmax,singletable);
TripletParameters & q = tripletParameters[i*nelements*nelements+k*nelements+j];
q.P = new TabularFunction;
(q.P)->set_values(nr,-cutmax,cutmax,singletable);
}
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
for (int i = 0; i < nelements; i++) {
TripletParameters & p = tripletParameters[i*nelements*nelements+j*nelements+k];
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
TripletParameters & q = tripletParameters[i*nelements*nelements+k*nelements+j];
q.P = new tabularFunction(nr,-cutmax,cutmax);
(q.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
}
}
}
}
if (eta == 3) {
@ -671,8 +711,8 @@ void PairPolymorphic::read_file(char *file)
reader->next_dvector(singletable, nr);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
p.P = new TabularFunction;
(p.P)->set_values(nr,-cutmax,cutmax,singletable);
}
}
delete[] singletable;
@ -683,8 +723,8 @@ void PairPolymorphic::read_file(char *file)
reader->next_dvector(singletable, ng);
}
MPI_Bcast(singletable,ng,MPI_DOUBLE,0,world);
p.G = new tabularFunction(ng,-1.0,1.0);
(p.G)->set_values(ng,-1.0,1.0,singletable,epsilon);
p.G = new TabularFunction;
(p.G)->set_values(ng,-1.0,1.0,singletable);
}
delete[] singletable;
singletable = new double[nx];
@ -694,8 +734,8 @@ void PairPolymorphic::read_file(char *file)
reader->next_dvector(singletable, nx);
}
MPI_Bcast(singletable,nx,MPI_DOUBLE,0,world);
p.F = new tabularFunction(nx,0.0,maxX);
(p.F)->set_values(nx,0.0,maxX,singletable,epsilon);
p.F = new TabularFunction;
(p.F)->set_values(nx,0.0,maxX,singletable);
}
delete[] singletable;
if (comm->me == 0) {
@ -739,28 +779,28 @@ void PairPolymorphic::setup_params()
n++;
}
for (i = 0; i < nelements-1; i++) {
for (j = i+1; j < nelements; j++) {
elem2param[match[i]][match[j]] = n;
elem2param[match[j]][match[i]] = n;
n++;
}
for (j = i+1; j < nelements; j++) {
elem2param[match[i]][match[j]] = n;
elem2param[match[j]][match[i]] = n;
n++;
}
}
// map atom triplet to parameter index
n = 0;
for (i = 0; i < nelements; i++)
for (j = 0; j < nelements; j++)
for (k = 0; k < nelements; k++) {
elem3param[match[i]][match[j]][match[k]] = n;
n++;
}
for (j = 0; j < nelements; j++)
for (k = 0; k < nelements; k++) {
elem3param[match[i]][match[j]][match[k]] = n;
n++;
}
// for debugging, call write_tables() to check the tabular functions
// if (comm->me == 0) {
// write_tables(51);
// }
// error->all(FLERR,"Test potential tables");
// for debugging, call write_tables() to check the tabular functions
#if defined(LMP_POLYMORPHIC_WRITE_TABLES)
if (comm->me == 0) write_tables(51);
error->all(FLERR,"Test potential tables");
#endif
}
/* ----------------------------------------------------------------------
@ -768,10 +808,10 @@ void PairPolymorphic::setup_params()
------------------------------------------------------------------------- */
void PairPolymorphic::attractive(PairParameters *p, PairParameters *q,
TripletParameters *trip,
double prefactor, double rij, double rik,
double *delrij, double *delrik,
double *fi, double *fj, double *fk)
TripletParameters *trip,
double prefactor, double rij, double rik,
double *delrij, double *delrik,
double *fi, double *fj, double *fk)
{
double rij_hat[3],rik_hat[3];
double rijinv,rikinv;
@ -788,20 +828,20 @@ void PairPolymorphic::attractive(PairParameters *p, PairParameters *q,
/* ---------------------------------------------------------------------- */
void PairPolymorphic::ters_zetaterm_d(double prefactor,
double *rij_hat, double rij,
double *rik_hat, double rik,
double *dri, double *drj, double *drk,
PairParameters *p, PairParameters *q,
TripletParameters *trip)
double *rij_hat, double rij,
double *rik_hat, double rik,
double *dri, double *drj, double *drk,
PairParameters *p, PairParameters *q,
TripletParameters *trip)
{
double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta;
double dcosdri[3],dcosdrj[3],dcosdrk[3];
cos_theta = dot3(rij_hat,rik_hat);
(q->W)->value(rik,fc,1,dfc,1);
(trip->P)->value(rij-(p->xi)*rik,ex_delr,1,ex_delr_d,1);
(trip->G)->value(cos_theta,gijk,1,gijk_d,1);
(q->W)->value2(rik,fc,1,dfc,1);
(trip->P)->value2(rij-(p->xi)*rik,ex_delr,1,ex_delr_d,1);
(trip->G)->value2(cos_theta,gijk,1,gijk_d,1);
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
@ -830,8 +870,8 @@ void PairPolymorphic::ters_zetaterm_d(double prefactor,
/* ---------------------------------------------------------------------- */
void PairPolymorphic::costheta_d(double *rij_hat, double rij,
double *rik_hat, double rik,
double *dri, double *drj, double *drk)
double *rik_hat, double rik,
double *dri, double *drj, double *drk)
{
// first element is devative wrt Ri, second wrt Rj, third wrt Rk
@ -846,107 +886,95 @@ void PairPolymorphic::costheta_d(double *rij_hat, double rij,
}
/* ---------------------------------------------------------------------- */
#if defined(LMP_POLYMORPHIC_WRITE_TABLES)
void PairPolymorphic::write_tables(int npts)
{
char tag[6] = "";
if (comm->me != 0) sprintf(tag,"%d",comm->me);
FILE* fp = nullptr;
double xmin,xmax,x,uf,vf,wf,pf,gf,ff,ufp,vfp,wfp,pfp,gfp,ffp;
char line[MAXLINE];
std::string filename;
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,"_UVW");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem2param[i][j];
PairParameters & pair = pairParameters[iparam_ij];
xmin = (pair.U)->get_xmin();
xmax = (pair.U)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.U)->value(x, uf, 1, ufp, 1);
(pair.V)->value(x, vf, 1, vfp, 1);
(pair.W)->value(x, wf, 1, wfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f %12.4f %12.4f %12.4f %12.4f \n",x,uf,vf,wf,ufp,vfp,wfp);
for (int j = 0; j < nelements; j++) {
filename = fmt::format("{}{}_UVW{}",elements[i],
elements[j],comm->me);
fp = fopen(filename.c_str(), "w");
int iparam_ij = elem2param[i][j];
PairParameters & pair = pairParameters[iparam_ij];
xmin = (pair.U)->get_xmin();
xmax = (pair.U)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.U)->value2(x, uf, 1, ufp, 1);
(pair.V)->value2(x, vf, 1, vfp, 1);
(pair.W)->value2(x, wf, 1, wfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f %12.4f %12.4f %12.4f %12.4f \n",x,uf,vf,wf,ufp,vfp,wfp);
}
fclose(fp);
}
fclose(fp);
}
}
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,elements[k]);
strcat(line,"_P");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem3param[i][j][k];
TripletParameters & pair = tripletParameters[iparam_ij];
xmin = (pair.P)->get_xmin();
xmax = (pair.P)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int n = 0; n < npts; n++) {
x = xmin + (xmax-xmin) * n / (npts-1);
(pair.P)->value(x, pf, 1, pfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,pf,pfp);
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) {
filename = fmt::format("{}{}{}_P{}",elements[i],elements[j],
elements[k],comm->me);
fp = fopen(filename.c_str(), "w");
int iparam_ij = elem3param[i][j][k];
TripletParameters & pair = tripletParameters[iparam_ij];
xmin = (pair.P)->get_xmin();
xmax = (pair.P)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int n = 0; n < npts; n++) {
x = xmin + (xmax-xmin) * n / (npts-1);
(pair.P)->value2(x, pf, 1, pfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,pf,pfp);
}
fclose(fp);
}
}
fclose(fp);
}
}
}
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,elements[k]);
strcat(line,"_G");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem3param[i][j][k];
TripletParameters & pair = tripletParameters[iparam_ij];
xmin = (pair.G)->get_xmin();
xmax = (pair.G)->get_xmax();
for (int n = 0; n < npts; n++) {
x = xmin + (xmax-xmin) * n / (npts-1);
(pair.G)->value(x, gf, 1, gfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,gf,gfp);
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) {
filename = fmt::format("{}{}{}_G{}",elements[i],elements[j],
elements[k],comm->me);
fp = fopen(filename.c_str(), "w");
int iparam_ij = elem3param[i][j][k];
TripletParameters & pair = tripletParameters[iparam_ij];
xmin = (pair.G)->get_xmin();
xmax = (pair.G)->get_xmax();
for (int n = 0; n < npts; n++) {
x = xmin + (xmax-xmin) * n / (npts-1);
(pair.G)->value2(x, gf, 1, gfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,gf,gfp);
}
fclose(fp);
}
}
fclose(fp);
}
}
}
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,"_F");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem2param[i][j];
PairParameters & pair = pairParameters[iparam_ij];
xmin = (pair.F)->get_xmin();
xmax = (pair.F)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.F)->value(x, ff, 1, ffp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,ff,ffp);
for (int j = 0; j < nelements; j++) {
filename = fmt::format("{}{}_F{}",elements[i],
elements[j],comm->me);
fp = fopen(filename.c_str(), "w");
int iparam_ij = elem2param[i][j];
PairParameters & pair = pairParameters[iparam_ij];
xmin = (pair.F)->get_xmin();
xmax = (pair.F)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.F)->value2(x, ff, 1, ffp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,ff,ffp);
}
fclose(fp);
}
fclose(fp);
}
}
}
#endif

263
src/MANYBODY/pair_polymorphic.h
View File

@ -24,10 +24,11 @@ PairStyle(polymorphic,PairPolymorphic)
#include <cmath>
namespace LAMMPS_NS {
// forward declaration
class TabularFunction;
class PairPolymorphic : public Pair {
public:
public:
PairPolymorphic(class LAMMPS *);
virtual ~PairPolymorphic();
@ -37,255 +38,25 @@ class PairPolymorphic : public Pair {
void init_style();
double init_one(int, int);
protected:
class tabularFunction {
public:
tabularFunction() {
size = 0;
xmin = 0.0;
xmax = 0.0;
xmaxsq = xmax*xmax;
vmax = 0.0;
xs = nullptr;
ys = nullptr;
ys1 = nullptr;
ys2 = nullptr;
ys3 = nullptr;
ys4 = nullptr;
ys5 = nullptr;
ys6 = nullptr;
}
tabularFunction(int n) {
size = n;
xmin = 0.0;
xmax = 0.0;
xmaxsq = xmax*xmax;
xs = new double[n];
ys = new double[n];
ys1 = new double[n];
ys2 = new double[n];
ys3 = new double[n];
ys4 = new double[n];
ys5 = new double[n];
ys6 = new double[n];
}
tabularFunction(int n, double x1, double x2) {
size = n;
xmin = x1;
xmax = x2;
xmaxsq = xmax*xmax;
xs = new double[n];
ys = new double[n];
ys1 = new double[n];
ys2 = new double[n];
ys3 = new double[n];
ys4 = new double[n];
ys5 = new double[n];
ys6 = new double[n];
}
virtual ~tabularFunction() {
delete [] xs;
delete [] ys;
delete [] ys1;
delete [] ys2;
delete [] ys3;
delete [] ys4;
delete [] ys5;
delete [] ys6;
}
void set_xrange(double x1, double x2) {
xmin = x1;
xmax = x2;
xmaxsq = xmax*xmax;
}
void set_values(int n, double x1, double x2, double * values, double epsilon)
{
int shrink = 1;
int ilo,ihi;
double vlo,vhi;
ilo = 0;
ihi = n-1;
for (int i = 0; i < n; i++) {
if (fabs(values[i]) <= epsilon) {
ilo = i;
} else {
break;
}
}
for (int i = n-1; i >= 0; i--) {
if (fabs(values[i]) <= epsilon) {
ihi = i;
} else {
break;
}
}
if (ihi < ilo) ihi = ilo;
vlo = values[ilo];
vhi = values[ilo];
for (int i = ilo; i <= ihi; i++) {
if (vlo > values[i]) vlo = values[i];
if (vhi < values[i]) vhi = values[i];
}
// shrink (remove near zero points) reduces cutoff radius, and therefore computational cost
// do not shrink when x2 < 1.1 (angular function) or x2 > 20.0 (non-radial function)
if (x2 < 1.1 || x2 > 20.0) {
shrink = 0;
}
// do not shrink when when list is abnormally small
if (ihi - ilo < 50) {
shrink = 0;
}
// shrink if it is a constant
if (vhi - vlo <= epsilon) {
// shrink = 1;
}
if (shrink == 0) {
ilo = 0;
ihi = n-1;
}
xmin = x1 + (x2-x1)/(n -1)*ilo;
xmax = xmin + (x2-x1)/(n -1)*(ihi-ilo);
xmaxsq = xmax*xmax;
n = ihi - ilo + 1;
resize(n);
for (int i = ilo; i <= ihi; i++) {
ys[i-ilo] = values[i];
}
initialize();
}
void value(double x, double &y, int ny, double &y1, int ny1)
{
double ps = (x - xmin) * rdx;
int ks = ps + 0.5;
if (ks > size-1) ks = size-1;
if (ks < 0 ) ks = 0;
ps = ps - ks;
if (ny) y = ((ys3[ks]*ps + ys2[ks])*ps + ys1[ks])*ps + ys[ks];
if (ny1) y1 = (ys6[ks]*ps + ys5[ks])*ps + ys4[ks];
}
void print_value()
{
printf("%d %f %f %f \n",size,xmin,xmax,rdx);
printf(" \n");
for (int i = 0; i < size; i++) {
printf("%f %f \n",xs[i],ys[i]);
}
}
double get_xmin() {
return xmin;
}
double get_xmax() {
return xmax;
}
double get_xmaxsq() {
return xmaxsq;
}
double get_rdx() {
return rdx;
}
double get_vmax() {
return vmax;
}
protected:
void resize(int n) {
if (n != size) {
size = n;
delete [] xs;
xs = new double[n];
delete [] ys;
ys = new double[n];
delete [] ys1;
ys1 = new double[n];
delete [] ys2;
ys2 = new double[n];
delete [] ys3;
ys3 = new double[n];
delete [] ys4;
ys4 = new double[n];
delete [] ys5;
ys5 = new double[n];
delete [] ys6;
ys6 = new double[n];
}
}
void initialize() {
int n = size;
rdx = (xmax-xmin)/(n-1.0);
vmax = 0.0;
for (int i = 0; i < n; i++) {
if (fabs(ys[i]) > vmax) vmax = fabs(ys[i]);
}
for (int i = 0; i < n; i++) {
xs[i] = xmin+i*rdx;
}
rdx = 1.0 / rdx;
ys1[0] = ys[1] - ys[0];
ys1[1] = 0.5 * (ys[2] - ys[0]);
ys1[n-2] = 0.5 * (ys[n-1] - ys[n-3]);
ys1[n-1] = ys[n-1] - ys[n-2];
for (int i = 2; i < n-2; i++) {
ys1[i]=((ys[i-2]-ys[i+2])+ 8.0*(ys[i+1]-ys[i-1]))/12.0;
}
for (int i = 0; i < n-1; i++) {
ys2[i]=3.0*(ys[i+1]-ys[i])-2.0*ys1[i]-ys1[i+1];
ys3[i]=ys1[i]+ys1[i+1]-2.0*(ys[i+1]-ys[i]);
}
ys2[n-1]=0.0;
ys3[n-1]=0.0;
for (int i = 0; i < n; i++) {
ys4[i]=ys1[i]*rdx;
ys5[i]=2.0*ys2[i]*rdx;
ys6[i]=3.0*ys3[i]*rdx;
}
}
int size;
double xmin,xmax,xmaxsq,rdx,vmax;
double *ys, *ys1, *ys2, *ys3, *ys4, *ys5, *ys6;
double *xs;
};
protected:
struct PairParameters {
double cut;
double cutsq;
double xi;
class tabularFunction *U;
class tabularFunction *V;
class tabularFunction *W;
class tabularFunction *F;
PairParameters() {
cut = 0.0;
cutsq = 0.0;
xi = 1.0;
U = nullptr;
V = nullptr;
W = nullptr;
F = nullptr;
};
~PairParameters() {
delete U;
delete V;
delete W;
delete F;
}
TabularFunction *U;
TabularFunction *V;
TabularFunction *W;
TabularFunction *F;
PairParameters();
~PairParameters();
};
struct TripletParameters {
class tabularFunction *P;
class tabularFunction *G;
TripletParameters() {
P = nullptr;
G = nullptr;
};
~TripletParameters() {
delete P;
delete G;
}
TabularFunction *P;
TabularFunction *G;
TripletParameters();
~TripletParameters();
};
double epsilon;
@ -311,8 +82,9 @@ class PairPolymorphic : public Pair {
virtual void read_file(char *);
void setup_params();
#if defined(LMP_POLYMORPHIC_WRITE_TABLES)
void write_tables(int);
#endif
void attractive(PairParameters *, PairParameters *, TripletParameters *,
double, double, double, double *, double *, double *,
double *, double *);
@ -323,7 +95,6 @@ class PairPolymorphic : public Pair {
void costheta_d(double *, double, double *, double,
double *, double *, double *);
};
}
#endif