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

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4326 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2010-06-18 21:54:42 +00:00
parent c71f0cf9f8
commit b3de0db9db
7 changed files with 109 additions and 127 deletions
Download Patch File Download Diff File Expand all files Collapse all files

38
src/MOLECULE/angle_table.cpp
View File

@ -189,8 +189,8 @@ void AngleTable::settings(int narg, char **arg)
else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE; else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
else error->all("Unknown table style in angle style table"); else error->all("Unknown table style in angle style table");
n = atoi(arg[1]); tablength = force->inumeric(arg[1]);
nm1 = n - 1; if (tablength < 2) error->all("Illegal number of angle table entries");
// delete old tables, since cannot just change settings // delete old tables, since cannot just change settings
@ -282,7 +282,7 @@ double AngleTable::equilibrium_angle(int i)
void AngleTable::write_restart(FILE *fp) void AngleTable::write_restart(FILE *fp)
{ {
fwrite(&tabstyle,sizeof(int),1,fp); fwrite(&tabstyle,sizeof(int),1,fp);
fwrite(&n,sizeof(int),1,fp); fwrite(&tablength,sizeof(int),1,fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -293,11 +293,10 @@ void AngleTable::read_restart(FILE *fp)
{ {
if (comm->me == 0) { if (comm->me == 0) {
fread(&tabstyle,sizeof(int),1,fp); fread(&tabstyle,sizeof(int),1,fp);
fread(&n,sizeof(int),1,fp); fread(&tablength,sizeof(int),1,fp);
} }
MPI_Bcast(&tabstyle,1,MPI_DOUBLE,0,world); MPI_Bcast(&tabstyle,1,MPI_DOUBLE,0,world);
MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(&tablength,1,MPI_INT,0,world);
nm1 = n - 1;
allocate(); allocate();
} }
@ -451,7 +450,8 @@ void AngleTable::compute_table(Table *tb)
{ {
// delta = table spacing in angle for N-1 bins // delta = table spacing in angle for N-1 bins
tb->delta = PI/ nm1; int tlm1 = tablength-1;
tb->delta = PI/ tlm1;
tb->invdelta = 1.0/tb->delta; tb->invdelta = 1.0/tb->delta;
tb->deltasq6 = tb->delta*tb->delta / 6.0; tb->deltasq6 = tb->delta*tb->delta / 6.0;
@ -460,31 +460,31 @@ void AngleTable::compute_table(Table *tb)
// de,df values = delta values of e,f // de,df values = delta values of e,f
// ang,e,f are N in length so de,df arrays can compute difference // ang,e,f are N in length so de,df arrays can compute difference
tb->ang = (double *) memory->smalloc(n*sizeof(double),"angle:ang"); tb->ang = (double *) memory->smalloc(tablength*sizeof(double),"angle:ang");
tb->e = (double *) memory->smalloc(n*sizeof(double),"angle:e"); tb->e = (double *) memory->smalloc(tablength*sizeof(double),"angle:e");
tb->de = (double *) memory->smalloc(nm1*sizeof(double),"angle:de"); tb->de = (double *) memory->smalloc(tlm1*sizeof(double),"angle:de");
tb->f = (double *) memory->smalloc(n*sizeof(double),"angle:f"); tb->f = (double *) memory->smalloc(tablength*sizeof(double),"angle:f");
tb->df = (double *) memory->smalloc(nm1*sizeof(double),"angle:df"); tb->df = (double *) memory->smalloc(tlm1*sizeof(double),"angle:df");
tb->e2 = (double *) memory->smalloc(n*sizeof(double),"angle:e2"); tb->e2 = (double *) memory->smalloc(tablength*sizeof(double),"angle:e2");
tb->f2 = (double *) memory->smalloc(n*sizeof(double),"angle:f2"); tb->f2 = (double *) memory->smalloc(tablength*sizeof(double),"angle:f2");
double a; double a;
for (int i = 0; i < n; i++) { for (int i = 0; i < tablength; i++) {
a = i*tb->delta; a = i*tb->delta;
tb->ang[i] = a; tb->ang[i] = a;
tb->e[i] = splint(tb->afile,tb->efile,tb->e2file,tb->ninput,a); tb->e[i] = splint(tb->afile,tb->efile,tb->e2file,tb->ninput,a);
tb->f[i] = splint(tb->afile,tb->ffile,tb->f2file,tb->ninput,a); tb->f[i] = splint(tb->afile,tb->ffile,tb->f2file,tb->ninput,a);
} }
for (int i = 0; i < nm1; i++) { for (int i = 0; i < tlm1; i++) {
tb->de[i] = tb->e[i+1] - tb->e[i]; tb->de[i] = tb->e[i+1] - tb->e[i];
tb->df[i] = tb->f[i+1] - tb->f[i]; tb->df[i] = tb->f[i+1] - tb->f[i];
} }
double ep0 = - tb->f[0]; double ep0 = - tb->f[0];
double epn = - tb->f[nm1]; double epn = - tb->f[tlm1];
spline(tb->ang,tb->e,n,ep0,epn,tb->e2); spline(tb->ang,tb->e,tablength,ep0,epn,tb->e2);
spline(tb->ang,tb->f,n,tb->fplo,tb->fphi,tb->f2); spline(tb->ang,tb->f,tablength,tb->fplo,tb->fphi,tb->f2);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

2
src/MOLECULE/angle_table.h
View File

@ -38,7 +38,7 @@ class AngleTable : public Angle {
double single(int, int, int, int); double single(int, int, int, int);
private: private:
int tabstyle,n,nm1; int tabstyle,tablength;
double *theta0; double *theta0;
struct Table { struct Table {

38
src/MOLECULE/bond_table.cpp
View File

@ -141,8 +141,8 @@ void BondTable::settings(int narg, char **arg)
else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE; else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
else error->all("Unknown table style in bond style table"); else error->all("Unknown table style in bond style table");
n = atoi(arg[1]); tablength = force->inumeric(arg[1]);
nm1 = n - 1; if (tablength < 2) error->all("Illegal number of bond table entries");
// delete old tables, since cannot just change settings // delete old tables, since cannot just change settings
@ -224,7 +224,7 @@ double BondTable::equilibrium_distance(int i)
void BondTable::write_restart(FILE *fp) void BondTable::write_restart(FILE *fp)
{ {
fwrite(&tabstyle,sizeof(int),1,fp); fwrite(&tabstyle,sizeof(int),1,fp);
fwrite(&n,sizeof(int),1,fp); fwrite(&tablength,sizeof(int),1,fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -235,11 +235,10 @@ void BondTable::read_restart(FILE *fp)
{ {
if (comm->me == 0) { if (comm->me == 0) {
fread(&tabstyle,sizeof(int),1,fp); fread(&tabstyle,sizeof(int),1,fp);
fread(&n,sizeof(int),1,fp); fread(&tablength,sizeof(int),1,fp);
} }
MPI_Bcast(&tabstyle,1,MPI_DOUBLE,0,world); MPI_Bcast(&tabstyle,1,MPI_DOUBLE,0,world);
MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(&tablength,1,MPI_INT,0,world);
nm1 = n - 1;
allocate(); allocate();
} }
@ -374,8 +373,9 @@ void BondTable::spline_table(Table *tb)
void BondTable::compute_table(Table *tb) void BondTable::compute_table(Table *tb)
{ {
// delta = table spacing for N-1 bins // delta = table spacing for N-1 bins
int tlm1 = tablength-1;
tb->delta = (tb->hi - tb->lo)/ nm1; tb->delta = (tb->hi - tb->lo)/ tlm1;
tb->invdelta = 1.0/tb->delta; tb->invdelta = 1.0/tb->delta;
tb->deltasq6 = tb->delta*tb->delta / 6.0; tb->deltasq6 = tb->delta*tb->delta / 6.0;
@ -384,31 +384,31 @@ void BondTable::compute_table(Table *tb)
// de,df values = delta values of e,f // de,df values = delta values of e,f
// r,e,f are N in length so de,df arrays can compute difference // r,e,f are N in length so de,df arrays can compute difference
tb->r = (double *) memory->smalloc(n*sizeof(double),"bond:r"); tb->r = (double *) memory->smalloc(tablength*sizeof(double),"bond:r");
tb->e = (double *) memory->smalloc(n*sizeof(double),"bond:e"); tb->e = (double *) memory->smalloc(tablength*sizeof(double),"bond:e");
tb->de = (double *) memory->smalloc(nm1*sizeof(double),"bond:de"); tb->de = (double *) memory->smalloc(tlm1*sizeof(double),"bond:de");
tb->f = (double *) memory->smalloc(n*sizeof(double),"bond:f"); tb->f = (double *) memory->smalloc(tablength*sizeof(double),"bond:f");
tb->df = (double *) memory->smalloc(nm1*sizeof(double),"bond:df"); tb->df = (double *) memory->smalloc(tlm1*sizeof(double),"bond:df");
tb->e2 = (double *) memory->smalloc(n*sizeof(double),"bond:e2"); tb->e2 = (double *) memory->smalloc(tablength*sizeof(double),"bond:e2");
tb->f2 = (double *) memory->smalloc(n*sizeof(double),"bond:f2"); tb->f2 = (double *) memory->smalloc(tablength*sizeof(double),"bond:f2");
double a; double a;
for (int i = 0; i < n; i++) { for (int i = 0; i < tablength; i++) {
a = tb->lo + i*tb->delta; a = tb->lo + i*tb->delta;
tb->r[i] = a; tb->r[i] = a;
tb->e[i] = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,a); tb->e[i] = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,a);
tb->f[i] = splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,a); tb->f[i] = splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,a);
} }
for (int i = 0; i < nm1; i++) { for (int i = 0; i < tlm1; i++) {
tb->de[i] = tb->e[i+1] - tb->e[i]; tb->de[i] = tb->e[i+1] - tb->e[i];
tb->df[i] = tb->f[i+1] - tb->f[i]; tb->df[i] = tb->f[i+1] - tb->f[i];
} }
double ep0 = - tb->f[0]; double ep0 = - tb->f[0];
double epn = - tb->f[nm1]; double epn = - tb->f[tlm1];
spline(tb->r,tb->e,n,ep0,epn,tb->e2); spline(tb->r,tb->e,tablength,ep0,epn,tb->e2);
spline(tb->r,tb->f,n,tb->fplo,tb->fphi,tb->f2); spline(tb->r,tb->f,tablength,tb->fplo,tb->fphi,tb->f2);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

2
src/MOLECULE/bond_table.h
View File

@ -38,7 +38,7 @@ class BondTable : public Bond {
double single(int, double, int, int); double single(int, double, int, int);
private: private:
int tabstyle,n,nm1; int tabstyle,tablength;
double *r0; double *r0;
struct Table { struct Table {

69
src/compute_dihedral_local.cpp
View File

@ -109,9 +109,8 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
{ {
int i,m,n,atom1,atom2,atom3,atom4; int i,m,n,atom1,atom2,atom3,atom4;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv;
double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; double s,c;
double c2mag,sin2,sc1,sc2,s1,s2,s12,c;
double *pbuf; double *pbuf;
double **x = atom->x; double **x = atom->x;
@ -148,71 +147,53 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
if (flag) { if (flag) {
// phi calculation from dihedral style OPLS // phi calculation from dihedral style harmonic
if (pflag >= 0) { if (pflag >= 0) {
vb1x = x[atom1][0] - x[atom2][0]; vb1x = x[atom1][0] - x[atom2][0];
vb1y = x[atom1][1] - x[atom2][1]; vb1y = x[atom1][1] - x[atom2][1];
vb1z = x[atom1][2] - x[atom2][2]; vb1z = x[atom1][2] - x[atom2][2];
domain->minimum_image(vb1x,vb1y,vb1z); domain->minimum_image(vb1x,vb1y,vb1z);
vb2x = x[atom3][0] - x[atom2][0]; vb2x = x[atom3][0] - x[atom2][0];
vb2y = x[atom3][1] - x[atom2][1]; vb2y = x[atom3][1] - x[atom2][1];
vb2z = x[atom3][2] - x[atom2][2]; vb2z = x[atom3][2] - x[atom2][2];
domain->minimum_image(vb2x,vb2y,vb2z); domain->minimum_image(vb2x,vb2y,vb2z);
vb2xm = -vb2x; vb2xm = -vb2x;
vb2ym = -vb2y; vb2ym = -vb2y;
vb2zm = -vb2z; vb2zm = -vb2z;
domain->minimum_image(vb2xm,vb2ym,vb2zm); domain->minimum_image(vb2xm,vb2ym,vb2zm);
vb3x = x[atom4][0] - x[atom3][0]; vb3x = x[atom4][0] - x[atom3][0];
vb3y = x[atom4][1] - x[atom3][1]; vb3y = x[atom4][1] - x[atom3][1];
vb3z = x[atom4][2] - x[atom3][2]; vb3z = x[atom4][2] - x[atom3][2];
domain->minimum_image(vb3x,vb3y,vb3z); domain->minimum_image(vb3x,vb3y,vb3z);
sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z);
sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z);
sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z);
rb1 = sqrt(sb1);
rb3 = sqrt(sb3);
c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3;
b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z;
b1mag = sqrt(b1mag2);
b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z;
b2mag = sqrt(b2mag2);
b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z;
b3mag = sqrt(b3mag2);
ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; ax = vb1y*vb2zm - vb1z*vb2ym;
r12c1 = 1.0 / (b1mag*b2mag); ay = vb1z*vb2xm - vb1x*vb2zm;
c1mag = ctmp * r12c1; az = vb1x*vb2ym - vb1y*vb2xm;
bx = vb3y*vb2zm - vb3z*vb2ym;
by = vb3z*vb2xm - vb3x*vb2zm;
bz = vb3x*vb2ym - vb3y*vb2xm;
ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; rasq = ax*ax + ay*ay + az*az;
r12c2 = 1.0 / (b2mag*b3mag); rbsq = bx*bx + by*by + bz*bz;
c2mag = ctmp * r12c2; rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
rg = sqrt(rgsq);
rginv = ra2inv = rb2inv = 0.0;
if (rg > 0) rginv = 1.0/rg;
if (rasq > 0) ra2inv = 1.0/rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq;
rabinv = sqrt(ra2inv*rb2inv);
sin2 = MAX(1.0 - c1mag*c1mag,0.0); c = (ax*bx + ay*by + az*bz)*rabinv;
sc1 = sqrt(sin2); s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
if (sc1 < SMALL) sc1 = SMALL;
sc1 = 1.0/sc1;
sin2 = MAX(1.0 - c2mag*c2mag,0.0);
sc2 = sqrt(sin2);
if (sc2 < SMALL) sc2 = SMALL;
sc2 = 1.0/sc2;
s1 = sc1 * sc1;
s2 = sc2 * sc2;
s12 = sc1 * sc2;
c = (c0 + c1mag*c2mag) * s12;
if (c > 1.0) c = 1.0; if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
pbuf[n] = 180.0*acos(c)/PI; pbuf[n] = 180.0*atan2(s,c)/PI;
} }
n += nvalues; n += nvalues;
} }

85
src/pair_table.cpp
View File

@ -76,6 +76,7 @@ void PairTable::compute(int eflag, int vflag)
Table *tb; Table *tb;
union_int_float_t rsq_lookup; union_int_float_t rsq_lookup;
int tlm1 = tablength - 1;
evdwl = 0.0; evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag); if (eflag || vflag) ev_setup(eflag,vflag);
@ -127,19 +128,19 @@ void PairTable::compute(int eflag, int vflag)
if (tabstyle == LOOKUP) { if (tabstyle == LOOKUP) {
itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta); itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta);
if (itable >= nm1) if (itable >= tlm1)
error->one("Pair distance > table outer cutoff"); error->one("Pair distance > table outer cutoff");
fpair = factor_lj * tb->f[itable]; fpair = factor_lj * tb->f[itable];
} else if (tabstyle == LINEAR) { } else if (tabstyle == LINEAR) {
itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta); itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta);
if (itable >= nm1) if (itable >= tlm1)
error->one("Pair distance > table outer cutoff"); error->one("Pair distance > table outer cutoff");
fraction = (rsq - tb->rsq[itable]) * tb->invdelta; fraction = (rsq - tb->rsq[itable]) * tb->invdelta;
value = tb->f[itable] + fraction*tb->df[itable]; value = tb->f[itable] + fraction*tb->df[itable];
fpair = factor_lj * value; fpair = factor_lj * value;
} else if (tabstyle == SPLINE) { } else if (tabstyle == SPLINE) {
itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta); itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta);
if (itable >= nm1) if (itable >= tlm1)
error->one("Pair distance > table outer cutoff"); error->one("Pair distance > table outer cutoff");
b = (rsq - tb->rsq[itable]) * tb->invdelta; b = (rsq - tb->rsq[itable]) * tb->invdelta;
a = 1.0 - b; a = 1.0 - b;
@ -220,8 +221,8 @@ void PairTable::settings(int narg, char **arg)
else if (strcmp(arg[0],"bitmap") == 0) tabstyle = BITMAP; else if (strcmp(arg[0],"bitmap") == 0) tabstyle = BITMAP;
else error->all("Unknown table style in pair_style command"); else error->all("Unknown table style in pair_style command");
n = force->inumeric(arg[1]); tablength = force->inumeric(arg[1]);
nm1 = n - 1; if (tablength < 2) error->all("Illegal number of pair table entries");
// delete old tables, since cannot just change settings // delete old tables, since cannot just change settings
@ -286,15 +287,13 @@ void PairTable::coeff(int narg, char **arg)
// match = 1 if don't need to spline read-in tables // match = 1 if don't need to spline read-in tables
// this is only the case if r values needed by final tables // this is only the case if r values needed by final tables
// exactly match r values read from file // exactly match r values read from file
// for tabstyle SPLINE, always need to build spline tables
tb->match = 0; tb->match = 0;
if (tabstyle == LINEAR && tb->ninput == n && if (tabstyle == LINEAR && tb->ninput == tablength &&
tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1; tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1;
if (tabstyle == SPLINE && tb->ninput == n && if (tabstyle == BITMAP && tb->ninput == 1 << tablength &&
tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1;
if (tabstyle == BITMAP && tb->ninput == 1 << n &&
tb->rflag == BMP && tb->rhi == tb->cut) tb->match = 1; tb->rflag == BMP && tb->rhi == tb->cut) tb->match = 1;
if (tb->rflag == BMP && tb->match == 0) if (tb->rflag == BMP && tb->match == 0)
error->all("Bitmapped table in file does not match requested table"); error->all("Bitmapped table in file does not match requested table");
@ -537,6 +536,8 @@ void PairTable::param_extract(Table *tb, char *line)
void PairTable::compute_table(Table *tb) void PairTable::compute_table(Table *tb)
{ {
int tlm1 = tablength-1;
// inner = inner table bound // inner = inner table bound
// cut = outer table bound // cut = outer table bound
// delta = table spacing in rsq for N-1 bins // delta = table spacing in rsq for N-1 bins
@ -545,7 +546,7 @@ void PairTable::compute_table(Table *tb)
if (tb->rflag) inner = tb->rlo; if (tb->rflag) inner = tb->rlo;
else inner = tb->rfile[0]; else inner = tb->rfile[0];
tb->innersq = inner*inner; tb->innersq = inner*inner;
tb->delta = (tb->cut*tb->cut - tb->innersq) / nm1; tb->delta = (tb->cut*tb->cut - tb->innersq) / tlm1;
tb->invdelta = 1.0/tb->delta; tb->invdelta = 1.0/tb->delta;
// direct lookup tables // direct lookup tables
@ -556,11 +557,11 @@ void PairTable::compute_table(Table *tb)
// e,f are never a match to read-in values, always computed via spline interp // e,f are never a match to read-in values, always computed via spline interp
if (tabstyle == LOOKUP) { if (tabstyle == LOOKUP) {
tb->e = (double *) memory->smalloc(nm1*sizeof(double),"pair:e"); tb->e = (double *) memory->smalloc(tlm1*sizeof(double),"pair:e");
tb->f = (double *) memory->smalloc(nm1*sizeof(double),"pair:f"); tb->f = (double *) memory->smalloc(tlm1*sizeof(double),"pair:f");
double r,rsq; double r,rsq;
for (int i = 0; i < nm1; i++) { for (int i = 0; i < tlm1; i++) {
rsq = tb->innersq + (i+0.5)*tb->delta; rsq = tb->innersq + (i+0.5)*tb->delta;
r = sqrt(rsq); r = sqrt(rsq);
tb->e[i] = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,r); tb->e[i] = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,r);
@ -577,14 +578,14 @@ void PairTable::compute_table(Table *tb)
// e,f can match read-in values, else compute via spline interp // e,f can match read-in values, else compute via spline interp
if (tabstyle == LINEAR) { if (tabstyle == LINEAR) {
tb->rsq = (double *) memory->smalloc(n*sizeof(double),"pair:rsq"); tb->rsq = (double *) memory->smalloc(tablength*sizeof(double),"pair:rsq");
tb->e = (double *) memory->smalloc(n*sizeof(double),"pair:e"); tb->e = (double *) memory->smalloc(tablength*sizeof(double),"pair:e");
tb->f = (double *) memory->smalloc(n*sizeof(double),"pair:f"); tb->f = (double *) memory->smalloc(tablength*sizeof(double),"pair:f");
tb->de = (double *) memory->smalloc(nm1*sizeof(double),"pair:de"); tb->de = (double *) memory->smalloc(tlm1*sizeof(double),"pair:de");
tb->df = (double *) memory->smalloc(nm1*sizeof(double),"pair:df"); tb->df = (double *) memory->smalloc(tlm1*sizeof(double),"pair:df");
double r,rsq; double r,rsq;
for (int i = 0; i < n; i++) { for (int i = 0; i < tablength; i++) {
rsq = tb->innersq + i*tb->delta; rsq = tb->innersq + i*tb->delta;
r = sqrt(rsq); r = sqrt(rsq);
tb->rsq[i] = rsq; tb->rsq[i] = rsq;
@ -597,7 +598,7 @@ void PairTable::compute_table(Table *tb)
} }
} }
for (int i = 0; i < nm1; i++) { for (int i = 0; i < tlm1; i++) {
tb->de[i] = tb->e[i+1] - tb->e[i]; tb->de[i] = tb->e[i+1] - tb->e[i];
tb->df[i] = tb->f[i+1] - tb->f[i]; tb->df[i] = tb->f[i+1] - tb->f[i];
} }
@ -612,16 +613,16 @@ void PairTable::compute_table(Table *tb)
// e,f can match read-in values, else compute via spline interp // e,f can match read-in values, else compute via spline interp
if (tabstyle == SPLINE) { if (tabstyle == SPLINE) {
tb->rsq = (double *) memory->smalloc(n*sizeof(double),"pair:rsq"); tb->rsq = (double *) memory->smalloc(tablength*sizeof(double),"pair:rsq");
tb->e = (double *) memory->smalloc(n*sizeof(double),"pair:e"); tb->e = (double *) memory->smalloc(tablength*sizeof(double),"pair:e");
tb->f = (double *) memory->smalloc(n*sizeof(double),"pair:f"); tb->f = (double *) memory->smalloc(tablength*sizeof(double),"pair:f");
tb->e2 = (double *) memory->smalloc(n*sizeof(double),"pair:e2"); tb->e2 = (double *) memory->smalloc(tablength*sizeof(double),"pair:e2");
tb->f2 = (double *) memory->smalloc(n*sizeof(double),"pair:f2"); tb->f2 = (double *) memory->smalloc(tablength*sizeof(double),"pair:f2");
tb->deltasq6 = tb->delta*tb->delta / 6.0; tb->deltasq6 = tb->delta*tb->delta / 6.0;
double r,rsq; double r,rsq;
for (int i = 0; i < n; i++) { for (int i = 0; i < tablength; i++) {
rsq = tb->innersq + i*tb->delta; rsq = tb->innersq + i*tb->delta;
r = sqrt(rsq); r = sqrt(rsq);
tb->rsq[i] = rsq; tb->rsq[i] = rsq;
@ -637,8 +638,8 @@ void PairTable::compute_table(Table *tb)
// ep0,epn = dE/dr at inner and at cut // ep0,epn = dE/dr at inner and at cut
double ep0 = - tb->f[0]; double ep0 = - tb->f[0];
double epn = - tb->f[nm1]; double epn = - tb->f[tlm1];
spline(tb->rsq,tb->e,n,ep0,epn,tb->e2); spline(tb->rsq,tb->e,tablength,ep0,epn,tb->e2);
// fp0,fpn = dh/dg at inner and at cut // fp0,fpn = dh/dg at inner and at cut
// h(r) = f(r)/r and g(r) = r^2 // h(r) = f(r)/r and g(r) = r^2
@ -657,17 +658,17 @@ void PairTable::compute_table(Table *tb)
} }
if (tb->fpflag && tb->cut == tb->rfile[tb->ninput-1]) fpn = if (tb->fpflag && tb->cut == tb->rfile[tb->ninput-1]) fpn =
(tb->fphi/tb->cut - tb->f[nm1]/(tb->cut*tb->cut)) / (2.0 * tb->cut); (tb->fphi/tb->cut - tb->f[tlm1]/(tb->cut*tb->cut)) / (2.0 * tb->cut);
else { else {
double rsq2 = tb->cut * tb->cut; double rsq2 = tb->cut * tb->cut;
double rsq1 = rsq2 - secant_factor*tb->delta; double rsq1 = rsq2 - secant_factor*tb->delta;
fpn = (tb->f[nm1] / sqrt(rsq2) - fpn = (tb->f[tlm1] / sqrt(rsq2) -
splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,sqrt(rsq1)) / splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,sqrt(rsq1)) /
sqrt(rsq1)) / (secant_factor*tb->delta); sqrt(rsq1)) / (secant_factor*tb->delta);
} }
for (int i = 0; i < n; i++) tb->f[i] /= sqrt(tb->rsq[i]); for (int i = 0; i < tablength; i++) tb->f[i] /= sqrt(tb->rsq[i]);
spline(tb->rsq,tb->f,n,fp0,fpn,tb->f2); spline(tb->rsq,tb->f,tablength,fp0,fpn,tb->f2);
} }
// bitmapped linear tables // bitmapped linear tables
@ -683,8 +684,8 @@ void PairTable::compute_table(Table *tb)
// linear lookup tables of length ntable = 2^n // linear lookup tables of length ntable = 2^n
// stored value = value at lower edge of bin // stored value = value at lower edge of bin
init_bitmap(inner,tb->cut,n,masklo,maskhi,tb->nmask,tb->nshiftbits); init_bitmap(inner,tb->cut,tablength,masklo,maskhi,tb->nmask,tb->nshiftbits);
int ntable = 1 << n; int ntable = 1 << tablength;
int ntablem1 = ntable - 1; int ntablem1 = ntable - 1;
tb->rsq = (double *) memory->smalloc(ntable*sizeof(double),"pair:rsq"); tb->rsq = (double *) memory->smalloc(ntable*sizeof(double),"pair:rsq");
@ -885,7 +886,7 @@ void PairTable::read_restart(FILE *fp)
void PairTable::write_restart_settings(FILE *fp) void PairTable::write_restart_settings(FILE *fp)
{ {
fwrite(&tabstyle,sizeof(int),1,fp); fwrite(&tabstyle,sizeof(int),1,fp);
fwrite(&n,sizeof(int),1,fp); fwrite(&tablength,sizeof(int),1,fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -896,11 +897,10 @@ void PairTable::read_restart_settings(FILE *fp)
{ {
if (comm->me == 0) { if (comm->me == 0) {
fread(&tabstyle,sizeof(int),1,fp); fread(&tabstyle,sizeof(int),1,fp);
fread(&n,sizeof(int),1,fp); fread(&tablength,sizeof(int),1,fp);
} }
MPI_Bcast(&tabstyle,1,MPI_DOUBLE,0,world); MPI_Bcast(&tabstyle,1,MPI_DOUBLE,0,world);
MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(&tablength,1,MPI_INT,0,world);
nm1 = n - 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -911,23 +911,24 @@ double PairTable::single(int i, int j, int itype, int jtype, double rsq,
{ {
int itable; int itable;
double fraction,value,a,b,phi; double fraction,value,a,b,phi;
int tlm1 = tablength - 1;
Table *tb = &tables[tabindex[itype][jtype]]; Table *tb = &tables[tabindex[itype][jtype]];
if (rsq < tb->innersq) error->one("Pair distance < table inner cutoff"); if (rsq < tb->innersq) error->one("Pair distance < table inner cutoff");
if (tabstyle == LOOKUP) { if (tabstyle == LOOKUP) {
itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta); itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta);
if (itable >= nm1) error->one("Pair distance > table outer cutoff"); if (itable >= tlm1) error->one("Pair distance > table outer cutoff");
fforce = factor_lj * tb->f[itable]; fforce = factor_lj * tb->f[itable];
} else if (tabstyle == LINEAR) { } else if (tabstyle == LINEAR) {
itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta); itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta);
if (itable >= nm1) error->one("Pair distance > table outer cutoff"); if (itable >= tlm1) error->one("Pair distance > table outer cutoff");
fraction = (rsq - tb->rsq[itable]) * tb->invdelta; fraction = (rsq - tb->rsq[itable]) * tb->invdelta;
value = tb->f[itable] + fraction*tb->df[itable]; value = tb->f[itable] + fraction*tb->df[itable];
fforce = factor_lj * value; fforce = factor_lj * value;
} else if (tabstyle == SPLINE) { } else if (tabstyle == SPLINE) {
itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta); itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta);
if (itable >= nm1) error->one("Pair distance > table outer cutoff"); if (itable >= tlm1) error->one("Pair distance > table outer cutoff");
b = (rsq - tb->rsq[itable]) * tb->invdelta; b = (rsq - tb->rsq[itable]) * tb->invdelta;
a = 1.0 - b; a = 1.0 - b;
value = a * tb->f[itable] + b * tb->f[itable+1] + value = a * tb->f[itable] + b * tb->f[itable+1] +

2
src/pair_table.h
View File

@ -40,7 +40,7 @@ class PairTable : public Pair {
void *extract(char *); void *extract(char *);
private: private:
int tabstyle,n,nm1; int tabstyle,tablength;
struct Table { struct Table {
int ninput,rflag,fpflag,match,ntablebits; int ninput,rflag,fpflag,match,ntablebits;
int nshiftbits,nmask; int nshiftbits,nmask;