jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Update Pair BOP

This commit is contained in:
Richard Berger 2020-06-04 14:36:25 -04:00
parent b2c4cce826
commit d9345a1652
No known key found for this signature in database
GPG Key ID: A9E83994E0BA0CAB
2 changed files with 341 additions and 413 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

752
src/MANYBODY/pair_bop.cpp
View File

@ -46,10 +46,11 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "tokenizer.h"
#include "potential_file_reader.h"
using namespace LAMMPS_NS;
#define MAXLINE 1024
#define EPSILON 1.0e-6
/* ---------------------------------------------------------------------- */
@ -294,7 +295,6 @@ void PairBOP::compute(int eflag, int vflag)
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int cnt1;
MPI_Comm_rank(world,&me);
inum = list->inum;
ilist = list->ilist;
@ -616,7 +616,6 @@ void PairBOP::coeff(int narg, char **arg)
{
int i,j;
int n = atom->ntypes;
MPI_Comm_rank(world,&me);
delete[] map;
map = new int[n+1];
@ -645,7 +644,7 @@ void PairBOP::coeff(int narg, char **arg)
read_table(arg[2]);
// match element names to BOP word types
if (me == 0) {
if (comm->me == 0) {
for (i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
@ -659,7 +658,7 @@ void PairBOP::coeff(int narg, char **arg)
}
MPI_Bcast(&map[1],atom->ntypes,MPI_INT,0,world);
if (me == 0) {
if (comm->me == 0) {
if (elements) {
for (i = 0; i < bop_types; i++) delete [] elements[i];
delete [] elements;
@ -1137,7 +1136,6 @@ double PairBOP::sigmaBo(int itmp, int jtmp)
nlocal = atom->nlocal;
ilist = list->ilist;
firstneigh = list->firstneigh;
MPI_Comm_rank(world,&me);
if(nb_sg==0) {
nb_sg=(maxneigh)*(maxneigh/2);
@ -4873,101 +4871,11 @@ double PairBOP::PiBo(int itmp, int jtmp)
}
/* ----------------------------------------------------------------------
read BOP potential file
allocate BOP tables
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void _noopt PairBOP::read_table(char *filename)
void PairBOP::allocate_tables()
{
int i,j,k,n,m;
int buf1,pass;
int nws,ws;
double buf2;
char s[MAXLINE],buf[2];
MPI_Comm_rank(world,&me);
if (me == 0) {
FILE *fp = force->open_potential(filename);
if (fp == NULL) {
char str[128];
snprintf(str,128,"Cannot open BOP potential file %s",filename);
error->one(FLERR,str);
}
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error); // skip first comment line
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%d",&bop_types);
elements = new char*[bop_types];
for(i=0;i<bop_types;i++) elements[i]=NULL;
for(i=0;i<bop_types;i++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
nws=0;
ws=1;
for(j=0;j<(int)strlen(s);j++) {
if(ws==1) {
if(isspace(s[j])) {
ws=1;
} else {
ws=0;
}
} else {
if(isspace(s[j])) {
ws=1;
nws++;
} else {
ws=0;
}
}
}
if(nws!=3){
error->all(FLERR,"Incorrect table format check for element types");
}
sscanf(s,"%d %lf %s",&buf1,&buf2,buf);
n= strlen(buf)+1;
elements[i] = new char[n];
strcpy(elements[i],buf);
}
nws=0;
ws=1;
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
for(j=0;j<(int)strlen(s);j++) {
if(ws==1) {
if(isspace(s[j])) {
ws=1;
} else {
ws=0;
}
} else {
if(isspace(s[j])) {
ws=1;
nws++;
} else {
ws=0;
}
}
}
if (nws==3) {
sscanf(s,"%d %d %d",&nr,&ntheta,&nBOt);
npower=2;
if(ntheta<=10) npower=ntheta;
} else if (nws==2) {
sscanf(s,"%d %d",&nr,&nBOt);
ntheta=0;
npower=3;
} else {
error->one(FLERR,"Unsupported BOP potential file format");
}
fclose(fp);
npairs=bop_types*(bop_types+1)/2;
}
MPI_Bcast(&nr,1,MPI_INT,0,world);
MPI_Bcast(&nBOt,1,MPI_INT,0,world);
MPI_Bcast(&ntheta,1,MPI_INT,0,world);
MPI_Bcast(&bop_types,1,MPI_INT,0,world);
MPI_Bcast(&npairs,1,MPI_INT,0,world);
MPI_Bcast(&npower,1,MPI_INT,0,world);
memory->destroy(pi_a);
memory->destroy(pro_delta);
memory->destroy(pi_delta);
@ -5011,330 +4919,350 @@ void _noopt PairBOP::read_table(char *filename)
memory->create(gfunc5,bop_types,bop_types,bop_types,ntheta,"BOP:gfunc5");
memory->create(gfunc6,bop_types,bop_types,bop_types,ntheta,"BOP:gfunc6");
memory->create(gpara,bop_types,bop_types,bop_types,npower+1,"BOP:gpara");
}
allocate();
if (me == 0) {
FILE *fp = force->open_potential(filename);
if (fp == NULL) {
char str[128];
snprintf(str,128,"Cannot open BOP potential file %s",filename);
error->one(FLERR,str);
}
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error); // skip first comment line
for(i=0;i<bop_types+2;i++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
}
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf%lf%lf",&small1,&small2,&small3g
,&small4,&small5,&small6,&small7);
for(i=0;i<bop_types;i++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf",&pi_p[i]);
}
cutmax=0;
for(i=0;i<npairs;i++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf",&rcut[i]);
if(rcut[i]>cutmax)
cutmax=rcut[i];
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf",&sigma_c[i],&sigma_a[i],&pi_c[i],&pi_a[i]);
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf",&sigma_delta[i],&pi_delta[i]);
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf",&sigma_f[i],&sigma_k[i],&small3[i]);
}
if(nws==3) {
for(i=0;i<bop_types;i++)
for(j=0;j<bop_types;j++)
for(k=j;k<bop_types;k++) {
if(npower<=2) {
for(m=0;m<ntheta;m++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&gfunc[j][i][k][n],&gfunc[j][i][k][n+1]
,&gfunc[j][i][k][n+2],&gfunc[j][i][k][n+3],&gfunc[j][i][k][n+4]);
n+=4;
/* ----------------------------------------------------------------------
read BOP potential file
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void _noopt PairBOP::read_table(char *filename)
{
if (comm->me == 0) {
try {
PotentialFileReader reader(lmp, filename, "BOP");
char * line = nullptr;
bop_types = reader.next_int();
elements = new char*[bop_types];
for(int i=0; i < bop_types; i++) {
line = reader.next_line(3);
ValueTokenizer values(line);
values.next_int();
values.next_double();
std::string name = values.next_string();
elements[i] = new char[name.length()+1];
strcpy(elements[i], name.c_str());
}
line = reader.next_line();
ValueTokenizer values(line);
int format = values.count();
switch(format) {
case 3:
nr = values.next_int();
ntheta = values.next_int();
nBOt = values.next_int();
if(ntheta <= 10)
npower = ntheta;
else
npower = 2;
break;
case 2:
nr = values.next_int();
nBOt = values.next_int();
ntheta = 0;
npower = 3;
break;
default:
error->one(FLERR,"Unsupported BOP potential file format");
}
npairs = bop_types*(bop_types+1)/2;
allocate_tables();
allocate();
line = reader.next_line(7);
values = ValueTokenizer(line);
small1 = values.next_double();
small2 = values.next_double();
small3g = values.next_double();
small4 = values.next_double();
small5 = values.next_double();
small6 = values.next_double();
small7 = values.next_double();
for(int i = 0; i < bop_types; i++) {
pi_p[i] = reader.next_double();
}
cutmax = 0.0;
for(int i = 0; i < npairs; i++) {
rcut[i] = reader.next_double();
if (rcut[i] > cutmax)
cutmax = rcut[i];
line = reader.next_line(4);
values = ValueTokenizer(line);
sigma_c[i] = values.next_double();
sigma_a[i] = values.next_double();
pi_c[i] = values.next_double();
pi_a[i] = values.next_double();
line = reader.next_line(2);
values = ValueTokenizer(line);
sigma_delta[i] = values.next_double();
pi_delta[i] = values.next_double();
line = reader.next_line(3);
values = ValueTokenizer(line);
sigma_f[i] = values.next_double();
sigma_k[i] = values.next_double();
small3[i] = values.next_double();
}
if(format == 3) {
for(int i = 0; i < bop_types; i++)
for(int j = 0; j < bop_types; j++)
for(int k = j; k < bop_types; k++) {
if(npower <= 2) {
reader.next_dvector(ntheta, &gfunc[j][i][k][0]);
} else {
reader.next_dvector(npower+1, &gpara[j][i][k][0]);
}
}
} else {
for(int i = 0; i < bop_types; i++)
for(int j = 0; j < bop_types; j++)
for(int k = 0; k < bop_types; k++) {
reader.next_dvector(3, &gpara[i][j][k][0]);
gpara[j][i][k][3] = 0;
}
}
for(int i = 0; i < npairs; i++) {
reader.next_dvector(nr, &pRepul[i][0]);
}
for(int i = 0; i < npairs; i++) {
reader.next_dvector(nr, &pBetaS[i][0]);
}
for(int i = 0; i < npairs; i++) {
reader.next_dvector(nr, &pBetaP[i][0]);
}
for(int i = 0; i < npairs; i++) {
reader.next_dvector(nBOt, &FsigBO[i][0]);
}
for(int i = 0; i < bop_types; i++) {
pro_delta[i] = reader.next_double();
}
for(int i = 0; i < bop_types; i++) {
pro[i] = reader.next_double();
}
for(int i=0;i<npairs;i++) {
rcut3[i]=0.0;
}
if(format == 3) {
for(int i = 0; i < npairs; i++) {
rcut3[i] = reader.next_double();
}
for(int i = 0; i < npairs; i++) {
reader.next_dvector(nr, &pLong[i][0]);
}
}
rcutall=0.0;
for(int i=0; i < npairs; i++) {
if(rcut[i]>rcutall)
rcutall=rcut[i];
if(rcut3[i]>rcutall)
rcutall=rcut3[i];
rcutsq[i]=rcut[i]*rcut[i];
rcutsq3[i]=rcut3[i]*rcut3[i];
dr[i]=rcut[i]/((double)nr-1.0);
rdr[i]=1.0/dr[i];
dr3[i]=rcut3[i]/((double)nr-1.0);
rdr3[i]=1.0/dr3[i];
}
rctroot=rcutall;
dtheta=2.0/((double)ntheta-1.0);
rdtheta=1.0/dtheta;
dBO=1.0/((double)nBOt-1.0);
rdBO=1.0/(double)dBO;
for (int i = 0; i < npairs; i++) {
pBetaS1[i][0]=pBetaS[i][1]-pBetaS[i][0];
pBetaS1[i][1]=0.5*(pBetaS[i][2]-pBetaS[i][0]);
pBetaS1[i][nr-2]=0.5*(pBetaS[i][nr-1]-pBetaS[i][nr-3]);
pBetaS1[i][nr-1]=pBetaS[i][nr-1]-pBetaS[i][nr-2];
pBetaP1[i][0]=pBetaP[i][1]-pBetaP[i][0];
pBetaP1[i][1]=0.5*(pBetaP[i][2]-pBetaP[i][0]);
pBetaP1[i][nr-2]=0.5*(pBetaP[i][nr-1]-pBetaP[i][nr-3]);
pBetaP1[i][nr-1]=pBetaP[i][nr-1]-pBetaP[i][nr-2];
pRepul1[i][0]=pRepul[i][1]-pRepul[i][0];
pRepul1[i][1]=0.5*(pRepul[i][2]-pRepul[i][0]);
pRepul1[i][nr-2]=0.5*(pRepul[i][nr-1]-pRepul[i][nr-3]);
pRepul1[i][nr-1]=pRepul[i][nr-1]-pRepul[i][nr-2];
FsigBO1[i][0]=FsigBO[i][1]-FsigBO[i][0];
FsigBO1[i][1]=0.5*(FsigBO[i][2]-FsigBO[i][0]);
FsigBO1[i][nBOt-2]=0.5*(FsigBO[i][nBOt-1]-FsigBO[i][nBOt-3]);
FsigBO1[i][nBOt-1]=FsigBO[i][nBOt-1]-FsigBO[i][nBOt-2];
pLong1[i][0]=pLong[i][1]-pLong[i][0];
pLong1[i][1]=0.5*(pLong[i][2]-pLong[i][0]);
pLong1[i][nBOt-2]=0.5*(pLong[i][nr-1]-pLong[i][nr-3]);
pLong1[i][nBOt-1]=pLong[i][nr-1]-pLong[i][nr-2];
for (int k = 2; k < nr-2; k++) {
pBetaS1[i][k]=((pBetaS[i][k-2]-pBetaS[i][k+2])
+8.0*(pBetaS[i][k+1]-pBetaS[i][k-1]))/12.0;
pBetaP1[i][k]=((pBetaP[i][k-2]-pBetaP[i][k+2])
+8.0*(pBetaP[i][k+1]-pBetaP[i][k-1]))/12.0;
pRepul1[i][k]=((pRepul[i][k-2]-pRepul[i][k+2])
+8.0*(pRepul[i][k+1]-pRepul[i][k-1]))/12.0;
pLong1[i][k]=((pLong[i][k-2]-pLong[i][k+2])
+8.0*(pLong[i][k+1]-pLong[i][k-1]))/12.0;
}
for (int k=2; k < nr-2; k++) {
FsigBO1[i][k]=((FsigBO[i][k-2]-FsigBO[i][k+2])
+8.0*(FsigBO[i][k+1]-FsigBO[i][k-1]))/12.0;
}
for (int k = 0; k < nr-1; k++) {
pBetaS2[i][k]=3.0*(pBetaS[i][k+1]-pBetaS[i][k])
-2.0*pBetaS1[i][k]-pBetaS1[i][k+1];
pBetaS3[i][k]=pBetaS1[i][k]+pBetaS1[i][k+1]
-2.0*(pBetaS[i][k+1]-pBetaS[i][k]);
pBetaP2[i][k]=3.0*(pBetaP[i][k+1]-pBetaP[i][k])
-2.0*pBetaP1[i][k]-pBetaP1[i][k+1];
pBetaP3[i][k]=pBetaP1[i][k]+pBetaP1[i][k+1]
-2.0*(pBetaP[i][k+1]-pBetaP[i][k]);
pRepul2[i][k]=3.0*(pRepul[i][k+1]-pRepul[i][k])
-2.0*pRepul1[i][k]-pRepul1[i][k+1];
pRepul3[i][k]=pRepul1[i][k]+pRepul1[i][k+1]
-2.0*(pRepul[i][k+1]-pRepul[i][k]);
pLong2[i][k]=3.0*(pLong[i][k+1]-pLong[i][k])
-2.0*pLong1[i][k]-pLong1[i][k+1];
pLong3[i][k]=pLong1[i][k]+pLong1[i][k+1]
-2.0*(pLong[i][k+1]-pLong[i][k]);
}
for (int k = 0; k < nBOt-1; k++) {
FsigBO2[i][k]=3.0*(FsigBO[i][k+1]-FsigBO[i][k])
-2.0*FsigBO1[i][k]-FsigBO1[i][k+1];
FsigBO3[i][k]=FsigBO1[i][k]+FsigBO1[i][k+1]
-2.0*(FsigBO[i][k+1]-FsigBO[i][k]);
}
pBetaS2[i][nr-1]=0.0;
pBetaS3[i][nr-1]=0.0;
pBetaP2[i][nr-1]=0.0;
pBetaP3[i][nr-1]=0.0;
pRepul2[i][nr-1]=0.0;
pRepul3[i][nr-1]=0.0;
pLong2[i][nr-1]=0.0;
pLong3[i][nr-1]=0.0;
FsigBO2[i][nBOt-1]=0.0;
FsigBO3[i][nBOt-1]=0.0;
for (int k=0; k < nr; k++) {
pBetaS4[i][k]=pBetaS1[i][k]/dr[i];
pBetaS5[i][k]=2.0*pBetaS2[i][k]/dr[i];
pBetaS6[i][k]=3.0*pBetaS3[i][k]/dr[i];
pBetaP4[i][k]=pBetaP1[i][k]/dr[i];
pBetaP5[i][k]=2.0*pBetaP2[i][k]/dr[i];
pBetaP6[i][k]=3.0*pBetaP3[i][k]/dr[i];
pRepul4[i][k]=pRepul1[i][k]/dr[i];
pRepul5[i][k]=2.0*pRepul2[i][k]/dr[i];
pRepul6[i][k]=3.0*pRepul3[i][k]/dr[i];
pLong4[i][k]=pLong1[i][k]/dr3[i];
pLong5[i][k]=2.0*pLong2[i][k]/dr3[i];
pLong6[i][k]=3.0*pLong3[i][k]/dr3[i];
}
for (int k=0; k < nBOt; k++) {
FsigBO4[i][k]=FsigBO1[i][k]/dBO;
FsigBO5[i][k]=2.0*FsigBO2[i][k]/dBO;
FsigBO6[i][k]=3.0*FsigBO3[i][k]/dBO;
}
}
if (npower <= 2) {
for (int i = 0; i < bop_types; i++) {
for (int j = 0; j < bop_types; j++) {
for (int k = j; k < bop_types; k++) {
gfunc1[j][i][k][0] = gfunc[j][i][k][1] - gfunc[j][i][k][0];
gfunc1[j][i][k][1] = 0.5 * (gfunc[j][i][k][2] - gfunc[j][i][k][0]);
gfunc1[j][i][k][ntheta - 2] = 0.5 * (gfunc[j][i][k][ntheta - 1] - gfunc[j][i][k][ntheta - 3]);
gfunc1[j][i][k][ntheta - 1] = 0.5 * (gfunc[j][i][k][ntheta - 1] - gfunc[j][i][k][ntheta - 2]);
for (int m = 2; m < ntheta - 2; m++) {
gfunc1[j][i][k][m] = ((gfunc[j][i][k][m - 2] - gfunc[j][i][k][m + 2]) +
8.0 * (gfunc[j][i][k][m + 1] - gfunc[j][i][k][m + 1] - gfunc[j][i][k][m - 1])) /
12.0;
}
for (int m = 0; m < ntheta - 1; m++) {
gfunc2[j][i][k][m] = 3.0 * (gfunc[j][i][k][m + 1] - gfunc[j][i][k][m]) -
2.0 * gfunc1[j][i][k][m] - gfunc1[j][i][k][m + 1];
gfunc3[j][i][k][m] = gfunc1[j][i][k][m] + gfunc1[j][i][k][m + 1] -
2.0 * (gfunc[j][i][k][m + 1] - gfunc[j][i][k][m]);
}
gfunc2[j][i][k][ntheta - 1] = 0.0;
gfunc3[j][i][k][ntheta - 1] = 0.0;
for (int m = 0; m < ntheta; m++) {
gfunc4[j][i][k][ntheta - 1] = gfunc1[j][i][k][m] / dtheta;
gfunc5[j][i][k][ntheta - 1] = 2.0 * gfunc2[j][i][k][m] / dtheta;
gfunc6[j][i][k][ntheta - 1] = 3.0 * gfunc3[j][i][k][m] / dtheta;
}
}
}
}
}
for (int i = 0; i < bop_types; i++) {
for (int j = 0; j < bop_types; j++) {
for (int k = 0; k < j; k++) {
if (npower <= 2) {
for (int n = 0; n < ntheta; n++) {
gfunc[j][i][k][n] = gfunc[k][i][j][n];
gfunc1[j][i][k][n] = gfunc1[k][i][j][n];
gfunc2[j][i][k][n] = gfunc2[k][i][j][n];
gfunc3[j][i][k][n] = gfunc3[k][i][j][n];
gfunc4[j][i][k][n] = gfunc4[k][i][j][n];
gfunc5[j][i][k][n] = gfunc5[k][i][j][n];
gfunc6[j][i][k][n] = gfunc6[k][i][j][n];
}
} else {
if(npower==3) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3]);
for (int n = 0; n < npower + 1; n++) {
gpara[j][i][k][n] = gpara[k][i][j][n];
}
else if(npower==4) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
}
else if(npower==5) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf",&gpara[j][i][k][5]);
}
else if(npower==6) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf",&gpara[j][i][k][5],&gpara[j][i][k][6]);
}
else if(npower==7) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf",&gpara[j][i][k][5],&gpara[j][i][k][6],&gpara[j][i][k][7]);
}
else if(npower==8) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf",&gpara[j][i][k][5],&gpara[j][i][k][6],&gpara[j][i][k][7],&gpara[j][i][k][8]);
}
else if(npower==9) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][5],&gpara[j][i][k][6],&gpara[j][i][k][7],&gpara[j][i][k][8],&gpara[j][i][k][9]);
}
else if(npower==10) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][5],&gpara[j][i][k][6],&gpara[j][i][k][7],&gpara[j][i][k][8],&gpara[j][i][k][9]);
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf",&gpara[j][i][k][10]);
}
}
}
} else {
for(i=0;i<bop_types;i++)
for(j=0;j<bop_types;j++)
for(k=0;k<bop_types;k++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf",&gpara[i][j][k][0],&gpara[i][j][k][1],&gpara[i][j][k][2]);
gpara[j][i][k][3]=0;
}
}
for(i=0;i<npairs;i++) {
for(j=0;j<nr;j++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&pRepul[i][j],&pRepul[i][j+1]
,&pRepul[i][j+2],&pRepul[i][j+3],&pRepul[i][j+4]);
j+=4;
}
}
for(i=0;i<npairs;i++) {
for(j=0;j<nr;j++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&pBetaS[i][j],&pBetaS[i][j+1]
,&pBetaS[i][j+2],&pBetaS[i][j+3],&pBetaS[i][j+4]);
j+=4;
}
}
for(i=0;i<npairs;i++) {
for(j=0;j<nr;j++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&pBetaP[i][j],&pBetaP[i][j+1]
,&pBetaP[i][j+2],&pBetaP[i][j+3],&pBetaP[i][j+4]);
j+=4;
}
}
for(i=0;i<npairs;i++) {
for(j=0;j<nBOt;j++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf%lf%lf%lf%lf",&FsigBO[i][j],&FsigBO[i][j+1]
,&FsigBO[i][j+2],&FsigBO[i][j+3],&FsigBO[i][j+4]);
j+=4;
}
}
for(i=0;i<bop_types;i++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf",&pro_delta[i]);
}
for(i=0;i<bop_types;i++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf",&pro[i]);
}
for(i=0;i<npairs;i++) {
rcut3[i]=0.0;
}
pass=0;
i=0;
if(nws==3) {
for(i=0;i<npairs;i++) {
utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
sscanf(s,"%lf",&rcut3[i]);
}
for(i=0;i<npairs;i++) {
for(j=0;j<nr;j++) {
pass=0;
while(fgets(s,MAXLINE,fp)!=NULL&&pass==0) {
sscanf(s,"%lf%lf%lf%lf%lf",&pLong[i][j],&pLong[i][j+1]
,&pLong[i][j+2],&pLong[i][j+3],&pLong[i][j+4]);
j+=4;
pass=1;
}
}
}
}
rcutall=0.0;
for(i=0;i<npairs;i++) {
if(rcut[i]>rcutall)
rcutall=rcut[i];
if(rcut3[i]>rcutall)
rcutall=rcut3[i];
rcutsq[i]=rcut[i]*rcut[i];
rcutsq3[i]=rcut3[i]*rcut3[i];
dr[i]=rcut[i]/((double)nr-1.0);
rdr[i]=1.0/dr[i];
dr3[i]=rcut3[i]/((double)nr-1.0);
rdr3[i]=1.0/dr3[i];
}
rctroot=rcutall;
dtheta=2.0/((double)ntheta-1.0);
rdtheta=1.0/dtheta;
dBO=1.0/((double)nBOt-1.0);
rdBO=1.0/(double)dBO;
for(i=0;i<npairs;i++) {
pBetaS1[i][0]=pBetaS[i][1]-pBetaS[i][0];
pBetaS1[i][1]=0.5*(pBetaS[i][2]-pBetaS[i][0]);
pBetaS1[i][nr-2]=0.5*(pBetaS[i][nr-1]-pBetaS[i][nr-3]);
pBetaS1[i][nr-1]=pBetaS[i][nr-1]-pBetaS[i][nr-2];
pBetaP1[i][0]=pBetaP[i][1]-pBetaP[i][0];
pBetaP1[i][1]=0.5*(pBetaP[i][2]-pBetaP[i][0]);
pBetaP1[i][nr-2]=0.5*(pBetaP[i][nr-1]-pBetaP[i][nr-3]);
pBetaP1[i][nr-1]=pBetaP[i][nr-1]-pBetaP[i][nr-2];
pRepul1[i][0]=pRepul[i][1]-pRepul[i][0];
pRepul1[i][1]=0.5*(pRepul[i][2]-pRepul[i][0]);
pRepul1[i][nr-2]=0.5*(pRepul[i][nr-1]-pRepul[i][nr-3]);
pRepul1[i][nr-1]=pRepul[i][nr-1]-pRepul[i][nr-2];
FsigBO1[i][0]=FsigBO[i][1]-FsigBO[i][0];
FsigBO1[i][1]=0.5*(FsigBO[i][2]-FsigBO[i][0]);
FsigBO1[i][nBOt-2]=0.5*(FsigBO[i][nBOt-1]-FsigBO[i][nBOt-3]);
FsigBO1[i][nBOt-1]=FsigBO[i][nBOt-1]-FsigBO[i][nBOt-2];
pLong1[i][0]=pLong[i][1]-pLong[i][0];
pLong1[i][1]=0.5*(pLong[i][2]-pLong[i][0]);
pLong1[i][nBOt-2]=0.5*(pLong[i][nr-1]-pLong[i][nr-3]);
pLong1[i][nBOt-1]=pLong[i][nr-1]-pLong[i][nr-2];
for(k=2;k<nr-2;k++) {
pBetaS1[i][k]=((pBetaS[i][k-2]-pBetaS[i][k+2])
+8.0*(pBetaS[i][k+1]-pBetaS[i][k-1]))/12.0;
pBetaP1[i][k]=((pBetaP[i][k-2]-pBetaP[i][k+2])
+8.0*(pBetaP[i][k+1]-pBetaP[i][k-1]))/12.0;
pRepul1[i][k]=((pRepul[i][k-2]-pRepul[i][k+2])
+8.0*(pRepul[i][k+1]-pRepul[i][k-1]))/12.0;
pLong1[i][k]=((pLong[i][k-2]-pLong[i][k+2])
+8.0*(pLong[i][k+1]-pLong[i][k-1]))/12.0;
}
for(k=2;k<nr-2;k++) {
FsigBO1[i][k]=((FsigBO[i][k-2]-FsigBO[i][k+2])
+8.0*(FsigBO[i][k+1]-FsigBO[i][k-1]))/12.0;
}
for(k=0;k<nr-1;k++) {
pBetaS2[i][k]=3.0*(pBetaS[i][k+1]-pBetaS[i][k])
-2.0*pBetaS1[i][k]-pBetaS1[i][k+1];
pBetaS3[i][k]=pBetaS1[i][k]+pBetaS1[i][k+1]
-2.0*(pBetaS[i][k+1]-pBetaS[i][k]);
pBetaP2[i][k]=3.0*(pBetaP[i][k+1]-pBetaP[i][k])
-2.0*pBetaP1[i][k]-pBetaP1[i][k+1];
pBetaP3[i][k]=pBetaP1[i][k]+pBetaP1[i][k+1]
-2.0*(pBetaP[i][k+1]-pBetaP[i][k]);
pRepul2[i][k]=3.0*(pRepul[i][k+1]-pRepul[i][k])
-2.0*pRepul1[i][k]-pRepul1[i][k+1];
pRepul3[i][k]=pRepul1[i][k]+pRepul1[i][k+1]
-2.0*(pRepul[i][k+1]-pRepul[i][k]);
pLong2[i][k]=3.0*(pLong[i][k+1]-pLong[i][k])
-2.0*pLong1[i][k]-pLong1[i][k+1];
pLong3[i][k]=pLong1[i][k]+pLong1[i][k+1]
-2.0*(pLong[i][k+1]-pLong[i][k]);
}
for(k=0;k<nBOt-1;k++) {
FsigBO2[i][k]=3.0*(FsigBO[i][k+1]-FsigBO[i][k])
-2.0*FsigBO1[i][k]-FsigBO1[i][k+1];
FsigBO3[i][k]=FsigBO1[i][k]+FsigBO1[i][k+1]
-2.0*(FsigBO[i][k+1]-FsigBO[i][k]);
}
pBetaS2[i][nr-1]=0.0;
pBetaS3[i][nr-1]=0.0;
pBetaP2[i][nr-1]=0.0;
pBetaP3[i][nr-1]=0.0;
pRepul2[i][nr-1]=0.0;
pRepul3[i][nr-1]=0.0;
pLong2[i][nr-1]=0.0;
pLong3[i][nr-1]=0.0;
FsigBO2[i][nBOt-1]=0.0;
FsigBO3[i][nBOt-1]=0.0;
for(k=0;k<nr;k++) {
pBetaS4[i][k]=pBetaS1[i][k]/dr[i];
pBetaS5[i][k]=2.0*pBetaS2[i][k]/dr[i];
pBetaS6[i][k]=3.0*pBetaS3[i][k]/dr[i];
pBetaP4[i][k]=pBetaP1[i][k]/dr[i];
pBetaP5[i][k]=2.0*pBetaP2[i][k]/dr[i];
pBetaP6[i][k]=3.0*pBetaP3[i][k]/dr[i];
pRepul4[i][k]=pRepul1[i][k]/dr[i];
pRepul5[i][k]=2.0*pRepul2[i][k]/dr[i];
pRepul6[i][k]=3.0*pRepul3[i][k]/dr[i];
pLong4[i][k]=pLong1[i][k]/dr3[i];
pLong5[i][k]=2.0*pLong2[i][k]/dr3[i];
pLong6[i][k]=3.0*pLong3[i][k]/dr3[i];
}
for(k=0;k<nBOt;k++) {
FsigBO4[i][k]=FsigBO1[i][k]/dBO;
FsigBO5[i][k]=2.0*FsigBO2[i][k]/dBO;
FsigBO6[i][k]=3.0*FsigBO3[i][k]/dBO;
}
}
if(npower<=2) {
for(i=0;i<bop_types;i++) {
for(j=0;j<bop_types;j++) {
for(k=j;k<bop_types;k++) {
gfunc1[j][i][k][0]=gfunc[j][i][k][1]-gfunc[j][i][k][0];
gfunc1[j][i][k][1]=0.5*(gfunc[j][i][k][2]-gfunc[j][i][k][0]);
gfunc1[j][i][k][ntheta-2]=0.5*(gfunc[j][i][k][ntheta-1]-gfunc[j][i][k][ntheta-3]);
gfunc1[j][i][k][ntheta-1]=0.5*(gfunc[j][i][k][ntheta-1]-gfunc[j][i][k][ntheta-2]);
for(m=2;m<ntheta-2;m++) {
gfunc1[j][i][k][m]=((gfunc[j][i][k][m-2]-gfunc[j][i][k][m+2])+
8.0*(gfunc[j][i][k][m+1]-gfunc[j][i][k][m+1]-gfunc[j][i][k][m-1]))/12.0;
}
for(m=0;m<ntheta-1;m++) {
gfunc2[j][i][k][m]=3.0*(gfunc[j][i][k][m+1]-gfunc[j][i][k][m])-
2.0*gfunc1[j][i][k][m]-gfunc1[j][i][k][m+1];
gfunc3[j][i][k][m]=gfunc1[j][i][k][m]+gfunc1[j][i][k][m+1]-
2.0*(gfunc[j][i][k][m+1]-gfunc[j][i][k][m]);
}
gfunc2[j][i][k][ntheta-1]=0.0;
gfunc3[j][i][k][ntheta-1]=0.0;
for(m=0;m<ntheta;m++) {
gfunc4[j][i][k][ntheta-1]=gfunc1[j][i][k][m]/dtheta;
gfunc5[j][i][k][ntheta-1]=2.0*gfunc2[j][i][k][m]/dtheta;
gfunc6[j][i][k][ntheta-1]=3.0*gfunc3[j][i][k][m]/dtheta;
}
}
}
}
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
} catch (FileReaderException & fre) {
error->one(FLERR, fre.what());
}
for(i=0;i<bop_types;i++) {
for(j=0;j<bop_types;j++) {
for(k=0;k<j;k++) {
if(npower<=2) {
for(n=0;n<ntheta;n++) {
gfunc[j][i][k][n]=gfunc[k][i][j][n];
gfunc1[j][i][k][n]=gfunc1[k][i][j][n];
gfunc2[j][i][k][n]=gfunc2[k][i][j][n];
gfunc3[j][i][k][n]=gfunc3[k][i][j][n];
gfunc4[j][i][k][n]=gfunc4[k][i][j][n];
gfunc5[j][i][k][n]=gfunc5[k][i][j][n];
gfunc6[j][i][k][n]=gfunc6[k][i][j][n];
}
} else {
for(n=0;n<npower+1;n++) {
gpara[j][i][k][n]=gpara[k][i][j][n];
}
}
}
}
}
fclose(fp);
}
MPI_Bcast(&nr,1,MPI_INT,0,world);
MPI_Bcast(&nBOt,1,MPI_INT,0,world);
MPI_Bcast(&ntheta,1,MPI_INT,0,world);
MPI_Bcast(&bop_types,1,MPI_INT,0,world);
MPI_Bcast(&npairs,1,MPI_INT,0,world);
MPI_Bcast(&npower,1,MPI_INT,0,world);
if (comm->me != 0){
allocate_tables();
allocate();
}
MPI_Bcast(&rdBO,1,MPI_DOUBLE,0,world);
MPI_Bcast(&dBO,1,MPI_DOUBLE,0,world);
MPI_Bcast(&rdtheta,1,MPI_DOUBLE,0,world);

2
src/MANYBODY/pair_bop.h
View File

@ -41,7 +41,6 @@ class PairBOP : public Pair {
double memory_usage();
private:
int me;
int maxneigh; // maximum size of neighbor list on this processor
int maxneigh3; // maximum size of neighbor list on this processor
int update_list; // check for changing maximum size of neighbor list
@ -207,6 +206,7 @@ class PairBOP : public Pair {
void read_table(char *);
void allocate();
void allocate_tables();
void create_pi(int);
void create_sigma(int);
void destroy_pi();