Homec4science
Diffusion lammps (lammm-master-deprecated) Browse

lammps/src/USER-MGPT/mgpt_readpot.cpp
38356949ba42lammm-master-deprecated

Last Change
Raw File

mgpt_readpot.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
This file is part of the MGPT implementation. See further comments
in pair_mgpt.cpp and pair_mgpt.h.
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "mgpt_splinetab.h"
#include "mgpt_readpot.h"
static double fgauss(double x,double al) {
return exp(-al * pow(x-1.0, 2));
}
static double hgauss(double x,double al) {
return (1.0 + al * pow(x-1.0, 2)) * exp(-al * pow(x-1.0, 2));
}
static double fl(double r,int mode,double rp,double p1,double al,double r0,double pn)
{
double term;
//double pn=1.0;
if (mode <= 4)
term = pow(rp/r, p1);
else
term = exp(-p1*(pow(r/rp, pn) - 1.0)/pn);
if (r <= r0) return term;
double quan = al*(r/r0 - 1.0)*(r/r0 - 1.0);
if(mode <= 2)
return term*exp(-quan);
else
return term*(1.0 + quan)*exp(-quan);
}
static int cmp_double(const void *ap,const void *bp) {
double a = *((const double *) ap);
double b = *((const double *) bp);
if(a < b)
return -1;
else if(a > b)
return 1;
else
return 0;
}
static void getparmindata(const char *potin_file,int nvol[1],double vol0[1],double x0[1],double x1[1]) {
int n,vsize;
double *volarr;
char metal[80],metalx[80];
int ipot,ipotx,mode,modex;
FILE *in = fopen(potin_file,"r");
char line[1024];
if(in == NULL) {
fprintf(stderr,"@%s:%d: Error reading potin file. Can not open file \'%s\'.\n",
__FILE__,__LINE__,potin_file);
exit(1);
}
vsize = 10;
volarr = (double *) malloc(sizeof(double) * vsize);
n = 0;
while(fgets(line,sizeof(line),in) != NULL) {
double zval,ivol,rws,mass;
double r0x,r1x,drx;
int nrx,i;
if(line[strspn(line," \t")] == '#') continue;
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);
}
} 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);
} 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);
}
}
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf %lf",&zval,&ivol,&rws,&mass);
if(n >= vsize) {
vsize = 2*vsize;
volarr = (double *) realloc(volarr,sizeof(double) * vsize);
}
volarr[n] = ivol;
n = n + 1;
for(i = 0; i<5; i++)
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf",&r0x,&r1x,&drx);
nrx = (int) ((r1x-r0x)/drx + 1.1); /* Really: 1+round((r1-r0)/dr) */
for(i = 0; i<nrx; i++)
fgets(line,sizeof(line),in);
}
fclose(in);
if(n == 0) {
fprintf(stderr,"@%s:%d: Invalid potin file \'%s\', no volume records.\n",
__FILE__,__LINE__,potin_file);
exit(1);
}
if(0) {
printf("Before sort:\n");
for(int i = 0; i<n; i++)
printf("%3d :: %.3f%s",i,volarr[i], (i%5==4) ? "\n" : " ");
printf("\n\n");
}
qsort(volarr,n,sizeof(double),cmp_double);
if(0) {
printf("After sort:\n");
for(int i = 0; i<n; i++)
printf("%3d :: %.3f%s",i,volarr[i], (i%5==4) ? "\n" : " ");
printf("\n\n");
}
nvol[0] = n;
vol0[0] = volarr[n/2];
x0[0] = pow(volarr[0]/vol0[0],1.0/3.0);
x1[0] = pow(volarr[n-1]/vol0[0],1.0/3.0);
free(volarr);
}
void potdata::readpot(const char *parmin_file,const char *potin_file,const double vol) {
FILE *in;
double x0,x1,dx;
int nx;
double r0x,r1x,drx;
int nrx;
char metalx[80];
int ipotx,modex; double pnx;
double vol0;
double *vatab,*vbtab,*vctab,*vdtab,*vetab,*p1tab,*altab,*vpairtab = 0;
double *r0rwstab,*evol0tab;
double (*C)[4];
double *y,*dy;
double x,dxdv;
double unused;
double zval,rws,ivol,r0rws,rcrws,rmrws;//,mass
int i,j;
int L;
char line[1024];
input_vol = vol;
/* Read potential data */
in = fopen(parmin_file,"r");
do {
fgets(line,sizeof(line),in);
} while(line[strspn(line," \t")] == '#');
/* Test to see whether this is a one-line or two-line version of parmin */
if(sscanf(line,"%lf %lf %lf %lf %d",&ddl[1],&ddl[2],&ddl[3],&ddl[4],&L) == 5) {
/* One-line version, call getparmindata to figure out volume table spacing. */
int nvol;
getparmindata(potin_file,&nvol,&vol0,&x0,&x1);
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);
}
} else {
/* Two-line version, reparse this line, and read second line */
sscanf(line,"%lf %lf %lf %lf",&x0,&x1,&dx,&vol0);
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf %lf %d",&ddl[1],&ddl[2],&ddl[3],&ddl[4],&L);
}
double rws_scale = pow(3.0*vol0/(16.0*atan(1.0)),1.0/3.0);
fclose(in);
lang = L+1;
lmax = 2*L+1;
double s = ddl[1],p = ddl[2],d = ddl[3],f = ddl[4];
double ss = s*s, pp = p*p, dd = d*d, ff = f*f;
anorm3 = s*ss + 2.0*( p*pp + d*dd + f*ff);
anorm4 = ss*ss + 2.0*(pp*pp + dd*dd + ff*ff);
/*
for(i = 1; i<=lmax; i++) {
for(j = 1; j<=lmax; j++)
del0.m[i][j] = 0.0;
del0[i][i] = 1.0;
}
Matrix::sz = lmax;
*/
nx = (int) ((x1-x0)/dx + 1.1); /* Really: 1+round((x1-x0)/dx) */
vatab = new double[nx];
vbtab = new double[nx];
vctab = new double[nx];
vdtab = new double[nx];
vetab = new double[nx];
p1tab = new double[nx];
altab = new double[nx];
r0rwstab = new double[nx];
evol0tab = new double[nx];
in = fopen(potin_file,"r");
int *tag = new int[nx];
for(i = 0; i<nx; i++) tag[i] = 0;
int ii;
for(ii = 0; ii<nx; ii++) {
do {
fgets(line,sizeof(line),in);
} while(line[strspn(line," \t")] == '#');
metalx[0] = 0;
/* 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);
}
if(modex <= 4) {
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);
}
} 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);
}
}
if(ii == 0) {
sscanf(line,"%s %d %d %lf",metal,&ipot,&mode,&pn);
if(modex <= 4) pn = pnx;
} 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);
}
}
//printf("LINE: %s\n",line);
//printf("metal = \'%s\' ipot = %d mode = %d\n",metalx,ipotx,modex);
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf %lf",&zval,&ivol,&rws,&mass);
/*{
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);
}*/
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);
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);
exit(1);
} else tag[i] = 1;
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf %lf",&r0rwstab[i],&altab[i],&rcrws,&rmrws);
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf",&p1tab[i],&unused,&evol0tab[i]);
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf %lf %lf",
&vatab[i],&vbtab[i],&vctab[i],&vdtab[i],&vetab[i]);
if(ipot == 1) {
vatab[i] *= vdtab[i];
vctab[i] *= vctab[i];
vetab[i] *= vetab[i];
}
fgets(line,sizeof(line),in);
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf",&r0x,&r1x,&drx);
nrx = (int) ((r1x-r0x)/drx + 1.1); /* Really: 1+round((r1-r0)/dr) */
if(ii == 0) {
r0 = r0x; r1 = r1x; nr = nrx;
vpairtab = new double[nx*nr];
} else {
/* Check that {r0,r1,dr,nr}x == {r0,r1,dr,nr} */
}
for(j = 0; j<nr; j++) {
double rj,ktan,dvdvol;
fgets(line,sizeof(line),in);
sscanf(line,"%lf %lf %lf %lf",
&rj,&vpairtab[i*nr+j],&ktan,&dvdvol);
{ /* Add screening and fl() part to pair energy table */
double al = altab[i];
double p1 = p1tab[i];
int bscreen = (al > 0.0);
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 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;
if(bscreen == 1 && rrws >= r0rws) {
double arg = rrws/r0rwstab[i];
double f;
if(mode <= 2) {
f = fgauss(arg,al);
}
else {
f = hgauss(arg,al);
}
fscr = f*f;
}
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 = NULL;
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);
}
}
}
}
fclose(in);
for(i = 0; i<nx; i++)
if(tag[i] == 0) {
printf("Volume point missing in table. i = %d\n",i);
exit(1);
}
/* Make table */
x = pow(vol/vol0,1.0/3.0);
dxdv = 1.0/(3.0*vol0*x*x);
C = new double[(nr > nx) ? nr : nx][4];
makespline(nx,1,vatab,C);
evalspline(nx-1,x0,x1,C,x,&va,&dva,&unused);
dva *= dxdv;
makespline(nx,1,vbtab,C);
evalspline(nx-1,x0,x1,C,x,&vb,&dvb,&unused);
dvb *= dxdv;
makespline(nx,1,vctab,C);
evalspline(nx-1,x0,x1,C,x,&vc,&dvc,&unused);
dvc *= dxdv;
makespline(nx,1,vdtab,C);
evalspline(nx-1,x0,x1,C,x,&vd,&dvd,&unused);
dvd *= dxdv;
makespline(nx,1,vetab,C);
evalspline(nx-1,x0,x1,C,x,&ve,&dve,&unused);
dve *= dxdv;
makespline(nx,1,p1tab,C);
evalspline(nx-1,x0,x1,C,x,&p1,&dp1,&unused);
dp1 *= dxdv;
makespline(nx,1,altab,C);
evalspline(nx-1,x0,x1,C,x,&al,&dal,&unused);
dal *= dxdv;
if(mode == 2 || mode == 4 || mode == 6) {
al = 125.0;
dal = 0.0;
}
{
double dr0rws;
makespline(nx,1,r0rwstab,C);
evalspline(nx-1,x0,x1,C,x,&r0rws,&dr0rws,&unused);
dr0rws *= dxdv;
rws = rws_scale*x;
r00 = r0rws * rws;
dr00 = dr0rws*rws + r0rws*rws_scale*dxdv;
rp = 1.8 * rws;
drp = 1.8 * rws_scale*dxdv;
}
makespline(nx,1,evol0tab,C);
evalspline(nx-1,x0,x1,C,x,&evol0,&devol0,&unused);
devol0 *= dxdv;
if(1) {
printf("%% READPOT PARAMETERS:\n");
printf("%% ddl = %15.5e %15.5e %15.5e %15.5e\n",ddl[1],ddl[2],ddl[3],ddl[4]);
printf("%% anorm3 = %15.5e anorm4 = %15.5e\n",anorm3,anorm4);
printf("%% x = %15.5e pn = %15.5e\n",x,pn);
printf("%% va = %15.5e dva = %15.5e\n",va,dva);
printf("%% vb = %15.5e dvb = %15.5e\n",vb,dvb);
printf("%% vc = %15.5e dvc = %15.5e\n",vc,dvc);
printf("%% vd = %15.5e dvd = %15.5e\n",vd,dvd);
printf("%% ve = %15.5e dve = %15.5e\n",ve,dve);
printf("%% p1 = %15.5e dp1 = %15.5e\n",p1,dp1);
printf("%% al = %15.5e dal = %15.5e\n",al,dal);
printf("%% rp = %15.5e drp = %15.5e\n",rp,drp);
printf("%% r00= %15.5e dr00= %15.5e\n",r00,dr00);
printf("\n");
}
y = new double[nr];
dy = new double[nr];
for(j = 0; j<nr; j++) {
double d2y;
makespline(nx,nr,&vpairtab[j],C);
evalspline(nx-1,x0,x1,C,x,&y[j],&dy[j],&d2y);
dy[j] *= dxdv;
}
vpair_spline = new double[nr-1][4];
dvpair_spline = new double[nr-1][4];
makespline(nr,1,y,vpair_spline);
makespline(nr,1,dy,dvpair_spline);
rcrit = rcrws * rws;
rmax = rmrws * rws;
delete[] dy;
delete[] y;
delete[] C;
delete[] evol0tab;
delete[] r0rwstab;
delete[] altab;
delete[] p1tab;
delete[] vetab;
delete[] vdtab;
delete[] vctab;
delete[] vbtab;
delete[] vatab;
}
/*
int main(int argc,char *argv[]) {
double vol = atof(argv[3]);
int n = 25,i;
printf("%% parmin = %s\n%% potin = %s\n%% vol = %15.5e\n",
argv[1],argv[2],vol);
readpot(argv[1],argv[2],vol);
for(i = 0; i<n; i++) {
double x,u,f,vir,dy,d2y;
x = r0 + i*(r1-r0)/(n-1);
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);
}
delete[] dvpair_spline;
delete[] vpair_spline;
return 0;
}
*/
c4science · Help