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

remove tabs

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-08-22 19:42:55 -04:00
parent 8b2e3ad3ee
commit 607e75f984
6 changed files with 1034 additions and 1034 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
src/MGPT/mgpt_linalg.h
View File

@ -41,13 +41,13 @@
#define const
#endif
typedef void (*trmul_fun) (const double * restrict A,
const double * restrict B,
double * restrict C);
const double * restrict B,
double * restrict C);
typedef void (*trtrace3_fun) (const double * restrict A,
const double * restrict B1,double * restrict t1,
const double * restrict B2,double * restrict t2,
const double * restrict B3,double * restrict t3);
const double * restrict B1,double * restrict t1,
const double * restrict B2,double * restrict t2,
const double * restrict B3,double * restrict t3);
#if defined(IBM_BG_SIMD) || defined(IBM_BGQ_SIMD)
#undef const
#endif

216
src/MGPT/mgpt_readpot.cpp
View File

@ -69,7 +69,7 @@ static void getparmindata(const char *potin_file,int nvol[1],double vol0[1],doub
if (in == nullptr) {
fprintf(stderr,"@%s:%d: Error reading potin file. Can not open file \'%s\'.\n",
__FILE__,__LINE__,potin_file);
__FILE__,__LINE__,potin_file);
exit(1);
}
@ -86,24 +86,24 @@ static void getparmindata(const char *potin_file,int nvol[1],double vol0[1],doub
if (n == 0) {
metal[0] = 0;
if (sscanf(line,"%s %d %d",metal,&ipot,&mode) != 3) {
fprintf(stderr,"@%s:%d: Error on potin file. line = %s\n",
__FILE__,__LINE__,line);
exit(1);
fprintf(stderr,"@%s:%d: Error on potin file. line = %s\n",
__FILE__,__LINE__,line);
exit(1);
}
} else {
metalx[0] = 0;
if (sscanf(line,"%s %d %d",metalx,&ipotx,&modex) != 3) {
fprintf(stderr,"@%s:%d: Error on potin file. line = %s\n",
__FILE__,__LINE__,line);
exit(1);
fprintf(stderr,"@%s:%d: Error on potin file. line = %s\n",
__FILE__,__LINE__,line);
exit(1);
} else if (strcmp(metal,metalx) != 0 || ipot != ipotx || mode != modex) {
fprintf(stderr,"@%s:%d: Error on potin file, parameter mismatch:\n"
" metal = \'%s\' ipot = %d mode = %d\n"
" metalx = \'%s\' ipotx = %d modex = %d\n",
__FILE__,__LINE__,
metal,ipot,mode,
metalx,ipotx,modex);
exit(1);
fprintf(stderr,"@%s:%d: Error on potin file, parameter mismatch:\n"
" metal = \'%s\' ipot = %d mode = %d\n"
" metalx = \'%s\' ipotx = %d modex = %d\n",
__FILE__,__LINE__,
metal,ipot,mode,
metalx,ipotx,modex);
exit(1);
}
}
@ -127,7 +127,7 @@ static void getparmindata(const char *potin_file,int nvol[1],double vol0[1],doub
if (n == 0) {
fprintf(stderr,"@%s:%d: Invalid potin file \'%s\', no volume records.\n",
__FILE__,__LINE__,potin_file);
__FILE__,__LINE__,potin_file);
exit(1);
}
@ -196,7 +196,7 @@ void potdata::readpot(const char *parmin_file,const char *potin_file,const doubl
dx = (x1-x0)/(nvol-1);
if (0) {
printf("getparmindata() ==> nvol = %d, vol0 = %.6f, x0= %.6f, x1 = %.6f, dx = %.6f\n",
nvol,vol0,x0,x1,dx);
nvol,vol0,x0,x1,dx);
}
} else {
/* Two-line version, reparse this line, and read second line */
@ -252,30 +252,30 @@ void potdata::readpot(const char *parmin_file,const char *potin_file,const doubl
/* Read element type, mode, and pn parameter */ {
int nf = sscanf(line,"%s %d %d %lf",metalx,&ipotx,&modex,&pnx);
if (nf < 3) {
printf("Error in %s() @ %s:%d: Inconsistency in potential input file (%s) "
"at record %d:\n"
" Expected at least three fields. Number of fields = %d\n",
__func__,__FILE__,__LINE__,potin_file,ii,
nf);
exit(1);
printf("Error in %s() @ %s:%d: Inconsistency in potential input file (%s) "
"at record %d:\n"
" Expected at least three fields. Number of fields = %d\n",
__func__,__FILE__,__LINE__,potin_file,ii,
nf);
exit(1);
}
if (modex <= 4) {
pnx = 1.0;
pnx = 1.0;
} else if (modex <= 6) {
if (nf != 4) {
printf("Error in %s() @ %s:%d: Inconsistency in potential input file (%s) "
"at record %d:\n"
" mode = %d, number of fields = %d\n",
__func__,__FILE__,__LINE__,potin_file,ii,
modex,nf);
exit(1);
}
if (nf != 4) {
printf("Error in %s() @ %s:%d: Inconsistency in potential input file (%s) "
"at record %d:\n"
" mode = %d, number of fields = %d\n",
__func__,__FILE__,__LINE__,potin_file,ii,
modex,nf);
exit(1);
}
} else {
printf("Error in %s() @ %s:%d: Inconsistency in potential input file (%s): "
"at record %d\n"
" Invalid mode. mode = %d\n",
__func__,__FILE__,__LINE__,potin_file,ii,
modex);
printf("Error in %s() @ %s:%d: Inconsistency in potential input file (%s): "
"at record %d\n"
" Invalid mode. mode = %d\n",
__func__,__FILE__,__LINE__,potin_file,ii,
modex);
}
}
@ -285,21 +285,21 @@ void potdata::readpot(const char *parmin_file,const char *potin_file,const doubl
} else {
/* Check that {metal,ipot,mode}x == {metal,ipot,mode} */
if (strcmp(metal,metalx) != 0 ||
ipotx != ipot ||
modex != mode ||
pnx != pn) {
printf("Error in %s() @ %s:%d: Inconsistency in potential input file (%s) "
"at record %d:\n"
"metalx != metal (%s != %s) or\n"
"ipotx != ipot (%d != %d) or\n"
"modex != mode (%d != %d) or\n"
"pnx != pn (%.3f != %.3f).\n",
__func__,__FILE__,__LINE__,potin_file,ii,
metalx,metal,
ipotx,ipot,
modex,mode,
pnx,pn);
exit(1);
ipotx != ipot ||
modex != mode ||
pnx != pn) {
printf("Error in %s() @ %s:%d: Inconsistency in potential input file (%s) "
"at record %d:\n"
"metalx != metal (%s != %s) or\n"
"ipotx != ipot (%d != %d) or\n"
"modex != mode (%d != %d) or\n"
"pnx != pn (%.3f != %.3f).\n",
__func__,__FILE__,__LINE__,potin_file,ii,
metalx,metal,
ipotx,ipot,
modex,mode,
pnx,pn);
exit(1);
}
}
//printf("LINE: %s\n",line);
@ -310,20 +310,20 @@ void potdata::readpot(const char *parmin_file,const char *potin_file,const doubl
double xi = x0 + i/((double) (nx-1)) * (x1-x0);
double volguess = vol0 * xi*xi*xi;
if (fabs(volguess/ivol - 1.0) > 1e-3)
printf("Wrong volume guess, i=%d volgues=%15.5e ivol=%15.5e\n",
i,volguess,ivol);
printf("Wrong volume guess, i=%d volgues=%15.5e ivol=%15.5e\n",
i,volguess,ivol);
}*/
double ifrac = (pow(ivol/vol0,1.0/3.0) - x0)/((x1-x0)/(nx-1));
i = (int) (ifrac + 0.1);
if (fabs(i - ifrac) > 0.01) {
printf("Volume point not in table... ii=%d i=%d ifrac=%15.5e vol=%15.5e\n",
ii,i,ifrac,ivol);
ii,i,ifrac,ivol);
printf("vol0 = %15.5e zval = %15.5e mass = %15.5e\n",vol0,zval,mass);
exit(1);
} else if (tag[i] == 1) {
printf("Duplicate volume point in table.... ii=%d i=%d ifrac=%15.5e vol=%15.5e\n",
ii,i,ifrac,ivol);
ii,i,ifrac,ivol);
exit(1);
} else tag[i] = 1;
@ -334,7 +334,7 @@ void potdata::readpot(const char *parmin_file,const char *potin_file,const doubl
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf %lf %lf",
&vatab[i],&vbtab[i],&vctab[i],&vdtab[i],&vetab[i]);
&vatab[i],&vbtab[i],&vctab[i],&vdtab[i],&vetab[i]);
if (ipot == 1) {
vatab[i] *= vdtab[i];
vctab[i] *= vctab[i];
@ -358,67 +358,67 @@ void potdata::readpot(const char *parmin_file,const char *potin_file,const doubl
double rj,ktan,dvdvol;
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf %lf",
&rj,&vpairtab[i*nr+j],&ktan,&dvdvol);
&rj,&vpairtab[i*nr+j],&ktan,&dvdvol);
{ /* Add screening and fl() part to pair energy table */
double al = altab[i];
double p1 = p1tab[i];
double al = altab[i];
double p1 = p1tab[i];
int bscreen = (al > 0.0);
int bscreen = (al > 0.0);
double xi = x0 + i/((double) (nx-1)) * (x1-x0);
double rws = rws_scale * xi;
double xi = x0 + i/((double) (nx-1)) * (x1-x0);
double rws = rws_scale * xi;
double r0rws = r0rwstab[i];
double r00 = r0rws*rws,rp = 1.8*rws;
if (bscreen == 0) r0rws = 10.0;
double alp = al,alm = al;
if (mode == 2 || mode == 4 || mode == 6) alm = 125.0;
al = alp;
double r0rws = r0rwstab[i];
double r00 = r0rws*rws,rp = 1.8*rws;
if (bscreen == 0) r0rws = 10.0;
double alp = al,alm = al;
if (mode == 2 || mode == 4 || mode == 6) alm = 125.0;
al = alp;
double r = r0 + j*(r1-r0)/(nr-1);
double r = r0 + j*(r1-r0)/(nr-1);
double rrws = r/rws;
//double rsqr = r*r;
// double fl(double r,int mode,double rp,double p1,double al,double r0)
double flr = fl(r,mode,rp,p1,al,r00,pn);
double fl2 = flr*flr;
double v2a = vatab[i]*fl2*fl2;
double v2b = vbtab[i]*fl2;
double fscr = 1.0;
double rrws = r/rws;
//double rsqr = r*r;
// double fl(double r,int mode,double rp,double p1,double al,double r0)
double flr = fl(r,mode,rp,p1,al,r00,pn);
double fl2 = flr*flr;
double v2a = vatab[i]*fl2*fl2;
double v2b = vbtab[i]*fl2;
double fscr = 1.0;
if (bscreen == 1 && rrws >= r0rws) {
double arg = rrws/r0rwstab[i];
double arg1 = arg - 1.0;
double arg12 = arg1*arg1;
double f,dp;
if (mode <= 2) {
f = fgauss(arg,al);
dp=2.*al*arg*arg1;
}
else {
f = hgauss(arg,al);
double arg13 = arg1*arg12;
dp=2.0*al*al*arg*arg13/(1.+al*arg12);
}
fscr = f*f;
}
if (bscreen == 1 && rrws >= r0rws) {
double arg = rrws/r0rwstab[i];
double arg1 = arg - 1.0;
double arg12 = arg1*arg1;
double f,dp;
if (mode <= 2) {
f = fgauss(arg,al);
dp=2.*al*arg*arg1;
}
else {
f = hgauss(arg,al);
double arg13 = arg1*arg12;
dp=2.0*al*al*arg*arg13/(1.+al*arg12);
}
fscr = f*f;
}
double vpair_tmp = vpairtab[i*nr+j];
vpairtab[i*nr+j] = vpairtab[i*nr+j]*fscr + v2a - v2b;
double vpair_tmp = vpairtab[i*nr+j];
vpairtab[i*nr+j] = vpairtab[i*nr+j]*fscr + v2a - v2b;
if (0) if (fabs(vol-ivol) < 0.01) {
static FILE *xfile = nullptr;
if (j == 0) {
xfile = fopen("mgpt5-pot.dat","w");
fprintf(xfile,"%%%% vol = %15.5e ivol = %15.5e i = %d ii = %d\n",
vol,ivol,i,ii);
}
fprintf(xfile,"%15.5e %15.5e %15.5e %15.5e %15.5e %20.10e\n",
r,vpair_tmp,fscr,v2a,v2b,flr);
if (j == nr-1) fclose(xfile);
}
if (0) if (fabs(vol-ivol) < 0.01) {
static FILE *xfile = nullptr;
if (j == 0) {
xfile = fopen("mgpt5-pot.dat","w");
fprintf(xfile,"%%%% vol = %15.5e ivol = %15.5e i = %d ii = %d\n",
vol,ivol,i,ii);
}
fprintf(xfile,"%15.5e %15.5e %15.5e %15.5e %15.5e %20.10e\n",
r,vpair_tmp,fscr,v2a,v2b,flr);
if (j == nr-1) fclose(xfile);
}
}
@ -544,7 +544,7 @@ int main(int argc,char *argv[]) {
int n = 25,i;
printf("%% parmin = %s\n%% potin = %s\n%% vol = %15.5e\n",
argv[1],argv[2],vol);
argv[1],argv[2],vol);
readpot(argv[1],argv[2],vol);
@ -555,7 +555,7 @@ int main(int argc,char *argv[]) {
evalspline(nr-1,r0,r1,vpair_spline,x,&u,&f,&d2y);
evalspline(nr-1,r0,r1,dvpair_spline,x,&vir,&dy,&d2y);
printf(" %15.5e %15.5e %15.5e %15.5e\n",
x,u,f,vir);
x,u,f,vir);
}

108
src/MGPT/mgpt_readpot.h
View File

@ -175,8 +175,8 @@ struct potdata2 {
p = strchr(s,'{');
if(p != nullptr) {
if(sscanf(p+1,"%d:%d:%d",i0,stride,i1) != 3) {
fprintf(stderr,"Error in template (\'%s\'), can not parse range.\n",nametemplate);
exit(1);
fprintf(stderr,"Error in template (\'%s\'), can not parse range.\n",nametemplate);
exit(1);
}
*p = '\0';
} else {
@ -205,13 +205,13 @@ struct potdata2 {
if(i0x != i0 || i1x != i1 || stridex != stride) {
fprintf(stderr,"Inconsistent templates. parmin_template=\'%s\', potin_template=\'%s\'\n",
parmin_template,potin_template);
parmin_template,potin_template);
exit(1);
}
if(i0 < 0 || i1 < i0 || stride <= 0 || (i1-i0)/stride+1 < 4) {
fprintf(stderr,"Improper temperature range. Need at least 4 temperature samples. "
"i0=%d,i1=%d,stride=%d,basename=\'%s\'\n",
i0,i1,stride,parmin_file);
"i0=%d,i1=%d,stride=%d,basename=\'%s\'\n",
i0,i1,stride,parmin_file);
exit(1);
}
@ -220,20 +220,20 @@ struct potdata2 {
if(parmin_suffix-1 == nullptr) {
fprintf(stderr,"No closing }. parmin_template=\'%s\'\n",
parmin_template);
parmin_template);
exit(1);
}
if(potin_suffix-1 == nullptr) {
fprintf(stderr,"No closing }. potin_template=\'%s\'\n",
potin_template);
potin_template);
exit(1);
}
printf("parmin_template = %s\n"
"parmin_file = %s\n"
"parmin_suffix = %s\n"
"T0=%d , T1=%d , stride=%d\n",
parmin_template,parmin_file,parmin_suffix,i0,i1,stride);
"parmin_file = %s\n"
"parmin_suffix = %s\n"
"T0=%d , T1=%d , stride=%d\n",
parmin_template,parmin_file,parmin_suffix,i0,i1,stride);
ntemp = (i1-i0)/stride + 1;
/*potdata **/potlist = new potdata[ntemp];
@ -245,27 +245,27 @@ struct potdata2 {
printf("Calling readpot(%s,%s,%.3f)\n",
parmin_file,potin_file,vol);
parmin_file,potin_file,vol);
potlist[k].readpot(parmin_file,potin_file,vol);
if(k > 0) {
if(potlist[k].nr != potlist[k-1].nr) {
fprintf(stderr,"nr differs between file %d and %d. Exiting.\n",
k,k-1);
exit(1);
}
if(potlist[k].nr != potlist[k-1].nr) {
fprintf(stderr,"nr differs between file %d and %d. Exiting.\n",
k,k-1);
exit(1);
}
if(potlist[k].r0 != potlist[k-1].r0) {
fprintf(stderr,"r0 differs between file %d and %d. Exiting.\n",
k,k-1);
exit(1);
}
if(potlist[k].r0 != potlist[k-1].r0) {
fprintf(stderr,"r0 differs between file %d and %d. Exiting.\n",
k,k-1);
exit(1);
}
if(potlist[k].r1 != potlist[k-1].r1) {
fprintf(stderr,"r1 differs between file %d and %d. Exiting.\n",
k,k-1);
exit(1);
}
if(potlist[k].r1 != potlist[k-1].r1) {
fprintf(stderr,"r1 differs between file %d and %d. Exiting.\n",
k,k-1);
exit(1);
}
}
}
tdeppot.r0 = potlist[0].r0;
@ -334,39 +334,39 @@ struct potdata2 {
tdeppot.dvpair = new double[sz][4][4];
/*
printf("vpair = %llx , dvpair = %llx",
(unsigned long long int) tdeppot.vpair,
(unsigned long long int) tdeppot.dvpair);
(unsigned long long int) tdeppot.vpair,
(unsigned long long int) tdeppot.dvpair);
printf(" @@@@@@@@@@@@@@ nr = %d\n",nr);
*/
for(int i = 0; i<nr-1; i++)
for(int j = 0; j<4; j++) {
/*
if(j == 5)
printf(" ############### i=%d\n",i);
*/
for(int j = 0; j<4; j++) {
/*
if(j == 5)
printf(" ############### i=%d\n",i);
*/
/* Make pair interaction interpolation functions */
for(int k = 0; k<ntemp; k++) {
if(0) if(i >= potlist[k].nr-1)
printf("Index error, local_nr=%d, k=%d, i=%d, nr=%d\n",nr,k,i,potlist[k].nr);
v[k] = potlist[k].vpair_spline[i][j];
}
makespline(ntemp,1,v,C);
/* Make pair interaction interpolation functions */
for(int k = 0; k<ntemp; k++) {
if(0) if(i >= potlist[k].nr-1)
printf("Index error, local_nr=%d, k=%d, i=%d, nr=%d\n",nr,k,i,potlist[k].nr);
v[k] = potlist[k].vpair_spline[i][j];
}
makespline(ntemp,1,v,C);
for(int k = 0; k<ntemp-1; k++)
for(int m = 0; m<4; m++)
tdeppot.vpair[k*(nr-1) + i][j][m] = C[k][m];
for(int k = 0; k<ntemp-1; k++)
for(int m = 0; m<4; m++)
tdeppot.vpair[k*(nr-1) + i][j][m] = C[k][m];
/* Make pair virial interpolation functions */
for(int k = 0; k<ntemp; k++)
v[k] = potlist[k].dvpair_spline[i][j];
makespline(ntemp,1,v,C);
/* Make pair virial interpolation functions */
for(int k = 0; k<ntemp; k++)
v[k] = potlist[k].dvpair_spline[i][j];
makespline(ntemp,1,v,C);
for(int k = 0; k<ntemp-1; k++)
for(int m = 0; m<4; m++)
tdeppot.dvpair[k*(nr-1) + i][j][m] = C[k][m];
for(int k = 0; k<ntemp-1; k++)
for(int m = 0; m<4; m++)
tdeppot.dvpair[k*(nr-1) + i][j][m] = C[k][m];
}
}
delete[] C;
delete[] v;
@ -391,7 +391,7 @@ struct potdata2 {
/*
printf("eval_pot nr=%d nt=%d\n",nr,nt);
printf("eval_pot r=%.3f T=%.3f k=%d i=%d\n",
r,T,k,i);
r,T,k,i);
printf("Tfrac=%.3f Tfrac-k=%.3f\n",Tfrac,Tfrac-k);
printf("rfrac=%.3f rfrac-i=%.3f\n",rfrac,rfrac-i);
*/
@ -401,7 +401,7 @@ struct potdata2 {
}
/*
printf("C coeff: %.3e %.3e %.3e %.3e\n",
C[0],C[1],C[2],C[3]);
C[0],C[1],C[2],C[3]);
*/
evalcubic(C,rfrac-i,e_p,f_p,&dd);
evalcubic(dC,rfrac-i,dedT_p,&dd,&dd);

2
src/MGPT/mgpt_splinetab.cpp
View File

@ -92,7 +92,7 @@ void evalcubic(double p[4],double x,double *y,double *dy,double *d2y) {
}
void evalspline(int n,double x0,double x1,double C[][4],
double x,double *y,double *dy,double *d2y) {
double x,double *y,double *dy,double *d2y) {
double xhat,t1,t2,t3;
double *p;
int idx;

2
src/MGPT/mgpt_splinetab.h
View File

@ -50,7 +50,7 @@ void makespline(int ntab,int stride,double tab[],double C[][4]);
The computational cost of this routine is O(1).
*/
void evalspline(int n,double x0,double x1,double C[][4],
double x,double *y,double *dy,double *d2y);
double x,double *y,double *dy,double *d2y);
/* Evaluate cubic polynomial represented by p in point x.

1730
src/MGPT/pair_mgpt.cpp
View File

File diff suppressed because it is too large Load Diff