tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2753 from akohlmey/update-pair-polymorphic

Browse Source 
Refactor pair style polymorphic to use the TabularFunction class shared with pair style bop
This commit is contained in:
Axel Kohlmeyer
2021-05-10 12:58:49 -04:00
committed by GitHub
parent 10d7a9f553 a4ea858bc2
commit 5895e960cb
10 changed files with 2315 additions and 2720 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
src/MANYBODY/pair_bop.cpp
View File

@ -42,6 +42,7 @@
#include "domain.h"
#include "error.h"
#include "force.h"
#include "math_special.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
@ -55,6 +56,9 @@
#include <string>
using namespace LAMMPS_NS;
using MathSpecial::square;
using MathSpecial::cube;
/* ---------------------------------------------------------------------- */
@ -188,7 +192,7 @@ PairBOP::~PairBOP()
void PairBOP::compute(int eflag, int vflag)
{
double minbox = MIN(MIN(domain->xprd, domain->yprd), domain->zprd);
if (minbox < 6.0*cutmax)
if (minbox-0.001 < 6.0*cutmax)
error->all(FLERR,"Pair style bop requires system dimension "
"of at least {:.4}",6.0*cutmax);
@ -308,9 +312,9 @@ void PairBOP::allocate()
tripletParameters = new TabularFunction[ntriples];
memory->create(elem2param,bop_types,bop_types,"BOP:elem2param");
memory->create(elem3param,bop_types,bop_types,bop_types,"BOP:elem3param");
bytes += npairs*sizeof(PairParameters) +
ntriples*sizeof(TabularFunction) + bop_types*bop_types*sizeof(int) +
bop_types*bop_types*bop_types*sizeof(int);
bytes += (double)npairs*sizeof(PairParameters) +
(double)ntriples*sizeof(TabularFunction) + square(bop_types)*sizeof(int) +
cube(bop_types)*sizeof(int);
memory->create(pi_a,npairs,"BOP:pi_a");
memory->create(pro_delta,bop_types,"BOP:pro_delta");
@ -403,8 +407,8 @@ void PairBOP::init_style()
// check that user sets comm->cutghostuser to 3x the max BOP cutoff
if (comm->cutghostuser < 3.0*cutmax)
error->all(FLERR,"Pair style bop requires a comm ghost cutoff "
if (comm->cutghostuser-0.001 < 3.0*cutmax)
error->all(FLERR,"Pair style bop requires setting a communication cutoff "
"of at least {:.4}",3.0*cutmax);
// need a full neighbor list and neighbors of ghosts

382
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)
@ -171,7 +211,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
if (eta == 1) {
iparam_ii = elem2param[itype][itype];
PairParameters & p = pairParameters[iparam_ii];
PairParameters &p = pairParameters[iparam_ii];
emb = (p.F)->get_vmax();
}
@ -191,13 +231,13 @@ void PairPolymorphic::compute(int eflag, int vflag)
r0 = sqrt(rsq);
iparam_ij = elem2param[itype][jtype];
PairParameters & p = pairParameters[iparam_ij];
PairParameters &p = pairParameters[iparam_ij];
// do not include the neighbor if get_vmax() <= epsilon because the function is near zero
if (eta == 1) {
if (emb > epsilon) {
iparam_jj = elem2param[jtype][jtype];
PairParameters & q = pairParameters[iparam_jj];
PairParameters &q = pairParameters[iparam_jj];
if (rsq < (q.W)->get_xmaxsq() && (q.W)->get_vmax() > epsilon) {
numneighW = numneighW + 1;
firstneighW[numneighW] = j;
@ -261,7 +301,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
if (emb > epsilon) {
iparam_ii = elem2param[itype][itype];
PairParameters & p = pairParameters[iparam_ii];
PairParameters &p = pairParameters[iparam_ii];
// accumulate bondorder zeta for each i-j interaction via loop over k
@ -272,7 +312,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
ktype = map[type[k]];
iparam_kk = elem2param[ktype][ktype];
PairParameters & q = pairParameters[iparam_kk];
PairParameters &q = pairParameters[iparam_kk];
(q.W)->value(drW[kk],wfac,1,fpair,0);
@ -299,7 +339,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
delr2[2] = -delzW[kk];
iparam_kk = elem2param[ktype][ktype];
PairParameters & q = pairParameters[iparam_kk];
PairParameters &q = pairParameters[iparam_kk];
(q.W)->value(drW[kk],wfac,0,fpair,1);
fpair = -prefactor*fpair/drW[kk];
@ -322,7 +362,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
jtype = map[type[j]];
iparam_ij = elem2param[itype][jtype];
PairParameters & p = pairParameters[iparam_ij];
PairParameters &p = pairParameters[iparam_ij];
delr1[0] = -delxV[jj];
delr1[1] = -delyV[jj];
@ -339,7 +379,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
if (j == k) continue;
ktype = map[type[k]];
iparam_ijk = elem3param[jtype][itype][ktype];
TripletParameters & trip = tripletParameters[iparam_ijk];
TripletParameters &trip = tripletParameters[iparam_ijk];
if ((trip.G)->get_vmax() <= epsilon) continue;
numneighW1 = numneighW1 + 1;
@ -354,7 +394,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
delr1[2]*delr2[2]) / (r1*r2);
iparam_ik = elem2param[itype][ktype];
PairParameters & q = pairParameters[iparam_ik];
PairParameters &q = pairParameters[iparam_ik];
(q.W)->value(r2,wfac,1,fpair,0);
(trip.P)->value(r1-(p.xi)*r2,pfac,1,fpair,0);
@ -388,7 +428,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
k = firstneighW[kk];
ktype = map[type[k]];
iparam_ijk = elem3param[jtype][itype][ktype];
TripletParameters & trip = tripletParameters[iparam_ijk];
TripletParameters &trip = tripletParameters[iparam_ijk];
delr2[0] = -delxW[kk];
delr2[1] = -delyW[kk];
@ -396,7 +436,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
r2 = drW[kk];
iparam_ik = elem2param[itype][ktype];
PairParameters & q = pairParameters[iparam_ik];
PairParameters &q = pairParameters[iparam_ik];
attractive(&p,&q,&trip,prefactor,r1,r2,delr1,delr2,fi,fj,fk);
@ -567,7 +607,7 @@ void PairPolymorphic::read_file(char *file)
tripletParameters = new TripletParameters[ntriple];
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
PairParameters &p = pairParameters[i];
values = reader->next_values(2);
p.cut = values.next_double();
p.cutsq = p.cut*p.cut;
@ -581,6 +621,7 @@ void PairPolymorphic::read_file(char *file)
MPI_Bcast(&nr, 1, MPI_INT, 0, world);
MPI_Bcast(&ng, 1, MPI_INT, 0, world);
MPI_Bcast(&nx, 1, MPI_INT, 0, world);
MPI_Bcast(&eta, 1, MPI_INT, 0, world);
MPI_Bcast(&maxX, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&npair, 1, MPI_INT, 0, world);
@ -601,110 +642,92 @@ void PairPolymorphic::read_file(char *file)
// start reading tabular functions
double * singletable = new double[nr];
for (int i = 0; i < npair; i++) { // U
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
PairParameters &p = pairParameters[i];
if (comm->me == 0) 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];
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
PairParameters &p = pairParameters[i];
if (comm->me == 0) 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];
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
PairParameters &p = pairParameters[i];
if (comm->me == 0) 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;
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
PairParameters &p = pairParameters[i];
if (p.cut > cutmax) cutmax = p.cut;
}
cutmaxsq = cutmax*cutmax;
if (eta != 3) {
for (int j = 0; j < nelements; j++) { // P
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
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+j];
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
TripletParameters &p = tripletParameters[i*nelements*nelements+j*nelements+j];
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) {
for (int i = 0; i < ntriple; i++) { // P
TripletParameters & p = tripletParameters[i];
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
TripletParameters &p = tripletParameters[i];
if (comm->me == 0) 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;
singletable = new double[ng];
for (int i = 0; i < ntriple; i++) { // G
TripletParameters & p = tripletParameters[i];
if (comm->me == 0) {
reader->next_dvector(singletable, ng);
}
TripletParameters &p = tripletParameters[i];
if (comm->me == 0) 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];
for (int i = 0; i < npair; i++) { // F
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
reader->next_dvector(singletable, nx);
}
PairParameters &p = pairParameters[i];
if (comm->me == 0) 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) {
delete reader;
}
if (comm->me == 0) delete reader;
// recalculate cutoffs of all params
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
PairParameters &p = pairParameters[i];
p.cut = (p.U)->get_xmax();
if (p.cut < (p.V)->get_xmax()) p.cut = (p.V)->get_xmax();
if (p.cut < (p.W)->get_xmax()) p.cut = (p.W)->get_xmax();
@ -714,7 +737,7 @@ void PairPolymorphic::read_file(char *file)
// set cutmax to max of all params
cutmax = 0.0;
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
PairParameters &p = pairParameters[i];
if (cutmax < p.cut) cutmax = p.cut;
}
cutmaxsq = cutmax*cutmax;
@ -739,28 +762,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 +791,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,11 +811,11 @@ 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];
@ -830,8 +853,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 +869,96 @@ 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)->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);
}
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)->value(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)->value(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)->value(x, ff, 1, ffp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,ff,ffp);
}
fclose(fp);
}
fclose(fp);
}
}
}
#endif

264
src/MANYBODY/pair_polymorphic.h
View File

@ -21,13 +21,13 @@ PairStyle(polymorphic,PairPolymorphic)
#define LMP_PAIR_POLYMORPHIC_H
#include "pair.h"
#include <cmath>
namespace LAMMPS_NS {
// forward declaration
class TabularFunction;
class PairPolymorphic : public Pair {
public:
public:
PairPolymorphic(class LAMMPS *);
virtual ~PairPolymorphic();
@ -37,255 +37,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 +81,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 +94,6 @@ class PairPolymorphic : public Pair {
void costheta_d(double *, double, double *, double,
double *, double *, double *);
};
}
#endif

223
src/tabular_function.cpp
View File

@ -18,71 +18,30 @@
using namespace LAMMPS_NS;
TabularFunction::TabularFunction()
: size(0), xmin(0.0), xmax(0.0), xmaxsq(0.0), rdx(0.0), vmax(0.0),
TabularFunction::TabularFunction() :
size(0), xmin(0.0), xmax(0.0), xmaxsq(0.0), rdx(0.0), vmax(0.0),
xs(nullptr), ys(nullptr), ys1(nullptr), ys2(nullptr), ys3(nullptr),
ys4(nullptr), ys5(nullptr), ys6(nullptr) {}
TabularFunction::TabularFunction(int n, double x1, double x2)
: TabularFunction()
ys4(nullptr), ys5(nullptr), ys6(nullptr)
{
set_xrange(x1, x2);
reset_size(n);
}
TabularFunction::~TabularFunction()
{
delete [] xs;
delete [] ys;
delete [] ys1;
delete [] ys2;
delete [] ys3;
delete [] ys4;
delete [] ys5;
delete [] ys6;
delete[] xs;
delete[] ys;
delete[] ys1;
delete[] ys2;
delete[] ys3;
delete[] ys4;
delete[] ys5;
delete[] ys6;
}
void TabularFunction::set_values(int n, double x1, double x2, double *values)
{
reset_size(n);
set_xrange(x1, x2);
memcpy(ys, values, n*sizeof(double));
initialize();
}
void TabularFunction::set_values(int n, double x1, double x2, double *values, double epsilon)
{
int ilo=0,ihi=n-1;
double vlo, vhi;
// get lo/hi boundaries where value < epsilon to shrink stored table
for (int i = ilo; ((fabs(values[i]) <= epsilon) && (i < n)); ++i)
ilo = i;
for (int i = ihi; ((fabs(values[i]) <= epsilon) && (i >= 0)); --i)
ihi = i;
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];
}
// do not shrink small tables
if (ihi - ilo < 50) {
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;
reset_size(n);
for (int i = ilo; i <= ihi; i++) {
ys[i-ilo] = values[i];
}
memcpy(ys, values, n * sizeof(double));
initialize();
}
@ -90,21 +49,21 @@ void TabularFunction::set_xrange(double x1, double x2)
{
xmin = x1;
xmax = x2;
xmaxsq = x2*x2;
xmaxsq = x2 * x2;
}
void TabularFunction::reset_size(int n)
{
if (n != size) {
size = n;
delete [] xs;
delete [] ys;
delete [] ys1;
delete [] ys2;
delete [] ys3;
delete [] ys4;
delete [] ys5;
delete [] ys6;
delete[] xs;
delete[] ys;
delete[] ys1;
delete[] ys2;
delete[] ys3;
delete[] ys4;
delete[] ys5;
delete[] ys6;
xs = new double[n];
ys = new double[n];
ys1 = new double[n];
@ -120,6 +79,9 @@ void TabularFunction::initialize()
{
int i;
rdx = (xmax - xmin) / (size - 1.0);
vmax = 0.0;
for (i = 0; i < size; i++)
if (fabs(ys[i]) > vmax) vmax = fabs(ys[i]);
for (i = 0; i < size; i++) xs[i] = xmin + i * rdx;
rdx = 1.0 / rdx;
ys1[0] = ys[1] - ys[0];
@ -140,140 +102,3 @@ void TabularFunction::initialize()
ys6[i] = 3.0 * ys3[i] * rdx;
}
}
#if 0
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]);
}
}
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;
};
#endif

24
src/tabular_function.h
View File

@ -18,11 +18,9 @@ namespace LAMMPS_NS {
class TabularFunction {
public:
TabularFunction();
TabularFunction(int, double, double);
virtual ~TabularFunction();
void set_values(int, double, double, double *);
void set_values(int, double, double, double *, double);
private:
int size;
@ -35,32 +33,18 @@ class TabularFunction {
public:
void value(double x, double &y, int ny, double &y1, int ny1)
{
double ps = (x - xmin) * rdx + 1.0;
int ks = ps;
if (ks > size - 1) ks = size - 1;
ps = ps - ks;
if (ps > 1.0) ps = 1.0;
if (ny)
y = ((ys3[ks - 1] * ps + ys2[ks - 1]) * ps + ys1[ks - 1]) * ps +
ys[ks - 1];
if (ny1) y1 = (ys6[ks - 1] * ps + ys5[ks - 1]) * ps + ys4[ks - 1];
}
void value2(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;
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];
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];
}
double get_xmin() const { return xmin; }
double get_xmax() const { return xmax; }
double get_xmaxsq() const { return xmaxsq; }
double get_rdx() const { return rdx; }
double get_vmax() { return vmax; }
};
} // namespace LAMMPS_NS

2
unittest/force-styles/tests/atomic-pair-polymorphic_eam.yaml
View File

@ -1,7 +1,7 @@
---
lammps_version: 10 Feb 2021
date_generated: Fri Feb 26 23:09:04 2021
epsilon: 1e-12
epsilon: 1e-14
prerequisites: ! |
pair polymorphic
pre_commands: ! |

2060
unittest/force-styles/tests/manybody-pair-bop.yaml
View File

File diff suppressed because it is too large Load Diff

2060
unittest/force-styles/tests/manybody-pair-bop_save.yaml
View File

File diff suppressed because it is too large Load Diff

2
unittest/force-styles/tests/manybody-pair-polymorphic_sw.yaml
View File

@ -1,7 +1,7 @@
---
lammps_version: 10 Feb 2021
date_generated: Fri Feb 26 23:09:17 2021
epsilon: 4e-12
epsilon: 1e-14
prerequisites: ! |
pair polymorphic
pre_commands: ! |

2
unittest/force-styles/tests/manybody-pair-polymorphic_tersoff.yaml
View File

@ -1,7 +1,7 @@
---
lammps_version: 10 Feb 2021
date_generated: Fri Feb 26 23:09:17 2021
epsilon: 1e-11
epsilon: 1e-14
prerequisites: ! |
pair polymorphic
pre_commands: ! |