Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/src/USER-MGPT/mgpt_readpot.h
354e20f43157lammm-devel

Last Change
Raw File

mgpt_readpot.h
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.
------------------------------------------------------------------------- */
#ifndef READPOT__
#define READPOT__
#include <cstdio>
#include "mgpt_splinetab.h"
struct potdata {
double va,vb,vc,vd,ve,p1,al,rp,r00,pn;
double dva,dvb,dvc,dvd,dve,dp1,dal,drp,dr00;
double evol0,devol0;
double (*vpair_spline)[4],(*dvpair_spline)[4];
double r0,r1;
int nr;
double mass,rcrit,rmax;
int lang,lmax;
double anorm3,anorm4;
double ddl[5];
int ipot,mode;
char metal[80];
double input_vol;
void readpot(const char *parmin_file,const char *potin_file,double vol);
void eval_pot(double r,double *e_p,double *f_p) {
double d2y;
evalspline(nr-1,r0,r1,vpair_spline,r,e_p,f_p,&d2y);
}
void eval_vir(double r,double *v_p) {
double dy,d2y;
evalspline(nr-1,r0,r1,dvpair_spline,r,v_p,&dy,&d2y);
}
};
struct potdata2 {
typedef double (*spline)[4];
spline va,vb,vc,vd,ve,p1,al,rp,r00;
spline dva,dvb,dvc,dvd,dve,dp1,dal,drp,dr00;
spline evol0,devol0;
double (*vpair)[4][4],(*dvpair)[4][4];
double r0,r1,T0,T1;
int nr,nt;
potdata *potlist;
double mass,rcrit,rmax;
int lang,lmax;
spline ddl[5];
int ipot,mode;
char metal[80];
double input_vol;
/* Functions to retrieve temperature dependent parameters */
double eval_tdep(spline y,double T) {
double f,df,d2f;
if(0) if(T != 3000.0)
printf("%s:%d: Error, T = %.3f\n",__FILE__,__LINE__,T);
evalspline(nt-1,T0,T1,y,T,&f,&df,&d2f);
return f;
}
double eval_tdepderiv(spline y,double T) {
double f,df,d2f;
if(0) if(T != 3000.0)
printf("%s:%d: Error, T = %.3f\n",__FILE__,__LINE__,T);
evalspline(nt-1,T0,T1,y,T,&f,&df,&d2f);
return df;
}
#define make_get(param) \
double get_##param(double T) { return eval_tdep(param,T); } \
double get_d##param(double T) { return eval_tdep(d##param,T); } \
double get_##param##_Tderiv(double T) { return eval_tdepderiv(param,T); }
/*
#define make_get(param) \
double get_##param(double T) { if(T != 3000.0) printf("%s:%d: Error, T = %.3f\n",__FILE__,__LINE__,T); return potlist[3].param; } \
double get_d##param(double T) { if(T != 3000.0) printf("%s:%d: Error, T = %.3f\n",__FILE__,__LINE__,T); return potlist[3].d##param; } \
double get_##param##_Tderiv(double T) { if(T != 3000.0) printf("%s:%d: Error, T = %.3f\n",__FILE__,__LINE__,T); return 0.0; }
*/
make_get(va) make_get(vb) make_get(vc) make_get(vd) make_get(ve)
make_get(p1) make_get(al) make_get(rp) make_get(r00) make_get(evol0)
#undef make_get
void get_anorm34(double T,double anorm3_p[1],double anorm4_p[1]) {
double
s = eval_tdep(ddl[1],T),
p = eval_tdep(ddl[2],T),
d = eval_tdep(ddl[3],T),
f = eval_tdep(ddl[4],T);
double ss = s*s, pp = p*p, dd = d*d, ff = f*f;
anorm3_p[0] = s*ss + 2.0*( p*pp + d*dd + f*ff);
anorm4_p[0] = ss*ss + 2.0*(pp*pp + dd*dd + ff*ff);
}
double get_anorm3(double T) {
double a3,a4;
get_anorm34(T,&a3,&a4);
return a3;
}
double get_anorm4(double T) {
double a3,a4;
get_anorm34(T,&a3,&a4);
return a4;
}
void get_anorm34_Tderiv(double T,double danorm3_p[1],double danorm4_p[1]) {
double
s = eval_tdep(ddl[1],T), ds = eval_tdepderiv(ddl[1],T),
p = eval_tdep(ddl[2],T), dp = eval_tdepderiv(ddl[2],T),
d = eval_tdep(ddl[3],T), d_d = eval_tdepderiv(ddl[3],T),
f = eval_tdep(ddl[4],T), df = eval_tdepderiv(ddl[4],T);
double ss = s*s, pp = p*p, dd = d*d, ff = f*f;
danorm3_p[0] = 3.0*ds*ss + 6.0*( dp*pp + d_d*dd + df*ff);
danorm4_p[0] = 4.0*ds*s*ss + 8.0*(dp*p*pp + d_d*d*dd + df*f*ff);
}
double get_anorm3_Tderiv(double T) {
double da3,da4;
get_anorm34_Tderiv(T,&da3,&da4);
return da3;
}
double get_anorm4_Tderiv(double T) {
double da3,da4;
get_anorm34_Tderiv(T,&da3,&da4);
return da4;
}
/* ... */
char * parsefname(const char *nametemplate,int *i0,int *i1,int *stride) {
char *s,*p;
if(0) {
s = new char[strlen(nametemplate)+1];
} else {
int len = 0;
while(nametemplate[len] != '\0') len = len + 1;
s = new char[len+1];
}
strcpy(s,nametemplate);
p = strchr(s,'{');
if(p != NULL) {
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);
}
*p = '\0';
} else {
*i0 = -1;
*i1 = -1;
*stride = 1;
}
return s;
}
spline maketempspline(int n,potdata data[],double *ptr) {
int stride = &(data[1].va) - &(data[0].va);
spline s = new double[n-1][4];
makespline(n,stride,ptr,s);
return s;
}
void readpot2(const char *parmin_template,const char *potin_template,double vol) {
int i0,i1,stride,i0x,i1x,stridex;
char *parmin_file = parsefname(parmin_template,&i0 ,&i1 ,&stride );
char *potin_file = parsefname( potin_template,&i0x,&i1x,&stridex);
int ntemp;
potdata2 &tdeppot = *this;
if(i0x != i0 || i1x != i1 || stridex != stride) {
fprintf(stderr,"Inconsistent templates. parmin_template=\'%s\', potin_template=\'%s\'\n",
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);
exit(1);
}
const char *parmin_suffix = strchr(parmin_template,'}')+1;
const char * potin_suffix = strchr( potin_template,'}')+1;
if(parmin_suffix-1 == NULL) {
fprintf(stderr,"No closing }. parmin_template=\'%s\'\n",
parmin_template);
exit(1);
}
if(potin_suffix-1 == NULL) {
fprintf(stderr,"No closing }. potin_template=\'%s\'\n",
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);
ntemp = (i1-i0)/stride + 1;
/*potdata **/potlist = new potdata[ntemp];
char *parend = parmin_file+strlen(parmin_file);
char *potend = potin_file +strlen( potin_file);
for(int k=0; k<ntemp; k++) {
sprintf(parend,"%d%s",i0+k*stride,parmin_suffix);
sprintf(potend,"%d%s",i0+k*stride,potin_suffix);
printf("Calling readpot(%s,%s,%.3f)\n",
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].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);
}
}
}
tdeppot.r0 = potlist[0].r0;
tdeppot.r1 = potlist[0].r1;
tdeppot.nr = potlist[0].nr;
tdeppot.T0 = i0;
tdeppot.T1 = i1;
tdeppot.nt = ntemp;
tdeppot.mass = potlist[0].mass;
tdeppot.rcrit = potlist[0].rcrit;
tdeppot.rmax = potlist[0].rmax;
tdeppot.lang = potlist[0].lang;
tdeppot.lmax = potlist[0].lmax;
tdeppot.ipot = potlist[0].ipot;
tdeppot.mode = potlist[0].mode;
tdeppot.input_vol = potlist[0].input_vol;
strncpy(tdeppot.metal,potlist[0].metal,sizeof(tdeppot.metal)/sizeof(char));
tdeppot.metal[sizeof(tdeppot.metal)/sizeof(char) - 1] = '\0';
delete[] parmin_file;
delete[] potin_file;
// Base parameters
tdeppot.va = maketempspline(ntemp,potlist,&(potlist[0].va));
tdeppot.vb = maketempspline(ntemp,potlist,&(potlist[0].vb));
tdeppot.vc = maketempspline(ntemp,potlist,&(potlist[0].vc));
tdeppot.vd = maketempspline(ntemp,potlist,&(potlist[0].vd));
tdeppot.ve = maketempspline(ntemp,potlist,&(potlist[0].ve));
tdeppot.p1 = maketempspline(ntemp,potlist,&(potlist[0].p1));
tdeppot.al = maketempspline(ntemp,potlist,&(potlist[0].al));
tdeppot.rp = maketempspline(ntemp,potlist,&(potlist[0].rp));
tdeppot.r00 = maketempspline(ntemp,potlist,&(potlist[0].r00));
tdeppot.evol0 = maketempspline(ntemp,potlist,&(potlist[0].evol0));
// Volume derivatives of base parameters
tdeppot.dva = maketempspline(ntemp,potlist,&(potlist[0].dva));
tdeppot.dvb = maketempspline(ntemp,potlist,&(potlist[0].dvb));
tdeppot.dvc = maketempspline(ntemp,potlist,&(potlist[0].dvc));
tdeppot.dvd = maketempspline(ntemp,potlist,&(potlist[0].dvd));
tdeppot.dve = maketempspline(ntemp,potlist,&(potlist[0].dve));
tdeppot.dp1 = maketempspline(ntemp,potlist,&(potlist[0].dp1));
tdeppot.dal = maketempspline(ntemp,potlist,&(potlist[0].dal));
tdeppot.drp = maketempspline(ntemp,potlist,&(potlist[0].drp));
tdeppot.dr00 = maketempspline(ntemp,potlist,&(potlist[0].dr00));
tdeppot.devol0 = maketempspline(ntemp,potlist,&(potlist[0].devol0));
tdeppot.ddl[0] = 0;
for(int k = 1; k<=4; k++)
tdeppot.ddl[k] = maketempspline(ntemp,potlist,&(potlist[0].ddl[k]));
{
double *v = new double[ntemp];
double (*C)[4] = new double[ntemp-1][4];
int sz = (nr-1)*(ntemp-1);
//printf("Allocation:: nr=%d ntemp=%d size=%d\n",nr,ntemp,sz);
tdeppot.vpair = new double[sz][4][4];
tdeppot.dvpair = new double[sz][4][4];
/*
printf("vpair = %llx , dvpair = %llx",
(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);
*/
/* 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];
/* 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];
}
delete[] C;
delete[] v;
}
}
void eval2Dspline(double fun[][4][4],double r,double T,double *e_p,double *f_p,double *dedT_p) {
double C[4],dC[4],dd;
double Tfrac = (T-T0)/(T1-T0) * (nt-1);
double rfrac = (r-r0)/(r1-r0) * (nr-1);
int k = (int) Tfrac;
int i = (int) rfrac;
int j;
if(k < 0) k = 0;
else if(k > nt-2) k = nt-2;
if(i < 0) i = 0;
else if(i > nr-2) i = nr-2;
/*
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);
printf("Tfrac=%.3f Tfrac-k=%.3f\n",Tfrac,Tfrac-k);
printf("rfrac=%.3f rfrac-i=%.3f\n",rfrac,rfrac-i);
*/
for(j = 0; j<4; j++) {
evalcubic(fun[k*(nr-1) + i][j],Tfrac-k,&C[j],&dC[j],&dd);
dC[j] = dC[j] * ((nt-1) / (T1-T0));
}
/*
printf("C coeff: %.3e %.3e %.3e %.3e\n",
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);
*f_p *= (nr-1) / (r1-r0);
}
void eval_pot(double r,double T,double *e_p,double *f_p,double *dedT_p) {
eval2Dspline(vpair,r,T,e_p,f_p,dedT_p);
}
void eval_vir(double r,double T,double *v_p) {
double vf,dvdT;
eval2Dspline(dvpair,r,T,v_p,&vf,&dvdT);
}
};
#endif
c4science · Help