From b3de0db9db300e7e3fc0e9b284e9b825603825af Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 21:54:42 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4326
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/MOLECULE/angle_table.cpp | 38 +++++++--------
src/MOLECULE/angle_table.h | 2 +-
src/MOLECULE/bond_table.cpp | 38 +++++++--------
src/MOLECULE/bond_table.h | 2 +-
src/compute_dihedral_local.cpp | 69 ++++++++++-----------------
src/pair_table.cpp | 85 +++++++++++++++++-----------------
src/pair_table.h | 2 +-
7 files changed, 109 insertions(+), 127 deletions(-)
diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp
index 1e6a5fcdfe..034db79e7e 100644
--- a/src/MOLECULE/angle_table.cpp
+++ b/src/MOLECULE/angle_table.cpp
@@ -189,8 +189,8 @@ void AngleTable::settings(int narg, char **arg)
else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
else error->all("Unknown table style in angle style table");
- n = atoi(arg[1]);
- nm1 = n - 1;
+ tablength = force->inumeric(arg[1]);
+ if (tablength < 2) error->all("Illegal number of angle table entries");
// delete old tables, since cannot just change settings
@@ -282,7 +282,7 @@ double AngleTable::equilibrium_angle(int i)
void AngleTable::write_restart(FILE *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) {
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(&n,1,MPI_INT,0,world);
- nm1 = n - 1;
+ MPI_Bcast(&tablength,1,MPI_INT,0,world);
allocate();
}
@@ -451,7 +450,8 @@ void AngleTable::compute_table(Table *tb)
{
// 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->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
// 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->e = (double *) memory->smalloc(n*sizeof(double),"angle:e");
- tb->de = (double *) memory->smalloc(nm1*sizeof(double),"angle:de");
- tb->f = (double *) memory->smalloc(n*sizeof(double),"angle:f");
- tb->df = (double *) memory->smalloc(nm1*sizeof(double),"angle:df");
- tb->e2 = (double *) memory->smalloc(n*sizeof(double),"angle:e2");
- tb->f2 = (double *) memory->smalloc(n*sizeof(double),"angle:f2");
+ tb->ang = (double *) memory->smalloc(tablength*sizeof(double),"angle:ang");
+ tb->e = (double *) memory->smalloc(tablength*sizeof(double),"angle:e");
+ tb->de = (double *) memory->smalloc(tlm1*sizeof(double),"angle:de");
+ tb->f = (double *) memory->smalloc(tablength*sizeof(double),"angle:f");
+ tb->df = (double *) memory->smalloc(tlm1*sizeof(double),"angle:df");
+ tb->e2 = (double *) memory->smalloc(tablength*sizeof(double),"angle:e2");
+ tb->f2 = (double *) memory->smalloc(tablength*sizeof(double),"angle:f2");
double a;
- for (int i = 0; i < n; i++) {
+ for (int i = 0; i < tablength; i++) {
a = i*tb->delta;
tb->ang[i] = 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);
}
- 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->df[i] = tb->f[i+1] - tb->f[i];
}
double ep0 = - tb->f[0];
- double epn = - tb->f[nm1];
- spline(tb->ang,tb->e,n,ep0,epn,tb->e2);
- spline(tb->ang,tb->f,n,tb->fplo,tb->fphi,tb->f2);
+ double epn = - tb->f[tlm1];
+ spline(tb->ang,tb->e,tablength,ep0,epn,tb->e2);
+ spline(tb->ang,tb->f,tablength,tb->fplo,tb->fphi,tb->f2);
}
/* ----------------------------------------------------------------------
diff --git a/src/MOLECULE/angle_table.h b/src/MOLECULE/angle_table.h
index a4895196d8..452ae13e13 100644
--- a/src/MOLECULE/angle_table.h
+++ b/src/MOLECULE/angle_table.h
@@ -38,7 +38,7 @@ class AngleTable : public Angle {
double single(int, int, int, int);
private:
- int tabstyle,n,nm1;
+ int tabstyle,tablength;
double *theta0;
struct Table {
diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp
index 7ffce66f1f..28e28f17f2 100644
--- a/src/MOLECULE/bond_table.cpp
+++ b/src/MOLECULE/bond_table.cpp
@@ -141,8 +141,8 @@ void BondTable::settings(int narg, char **arg)
else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
else error->all("Unknown table style in bond style table");
- n = atoi(arg[1]);
- nm1 = n - 1;
+ tablength = force->inumeric(arg[1]);
+ if (tablength < 2) error->all("Illegal number of bond table entries");
// delete old tables, since cannot just change settings
@@ -224,7 +224,7 @@ double BondTable::equilibrium_distance(int i)
void BondTable::write_restart(FILE *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) {
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(&n,1,MPI_INT,0,world);
- nm1 = n - 1;
+ MPI_Bcast(&tablength,1,MPI_INT,0,world);
allocate();
}
@@ -374,8 +373,9 @@ void BondTable::spline_table(Table *tb)
void BondTable::compute_table(Table *tb)
{
// 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->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
// 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->e = (double *) memory->smalloc(n*sizeof(double),"bond:e");
- tb->de = (double *) memory->smalloc(nm1*sizeof(double),"bond:de");
- tb->f = (double *) memory->smalloc(n*sizeof(double),"bond:f");
- tb->df = (double *) memory->smalloc(nm1*sizeof(double),"bond:df");
- tb->e2 = (double *) memory->smalloc(n*sizeof(double),"bond:e2");
- tb->f2 = (double *) memory->smalloc(n*sizeof(double),"bond:f2");
+ tb->r = (double *) memory->smalloc(tablength*sizeof(double),"bond:r");
+ tb->e = (double *) memory->smalloc(tablength*sizeof(double),"bond:e");
+ tb->de = (double *) memory->smalloc(tlm1*sizeof(double),"bond:de");
+ tb->f = (double *) memory->smalloc(tablength*sizeof(double),"bond:f");
+ tb->df = (double *) memory->smalloc(tlm1*sizeof(double),"bond:df");
+ tb->e2 = (double *) memory->smalloc(tablength*sizeof(double),"bond:e2");
+ tb->f2 = (double *) memory->smalloc(tablength*sizeof(double),"bond:f2");
double a;
- for (int i = 0; i < n; i++) {
+ for (int i = 0; i < tablength; i++) {
a = tb->lo + i*tb->delta;
tb->r[i] = 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);
}
- 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->df[i] = tb->f[i+1] - tb->f[i];
}
double ep0 = - tb->f[0];
- double epn = - tb->f[nm1];
- spline(tb->r,tb->e,n,ep0,epn,tb->e2);
- spline(tb->r,tb->f,n,tb->fplo,tb->fphi,tb->f2);
+ double epn = - tb->f[tlm1];
+ spline(tb->r,tb->e,tablength,ep0,epn,tb->e2);
+ spline(tb->r,tb->f,tablength,tb->fplo,tb->fphi,tb->f2);
}
/* ----------------------------------------------------------------------
diff --git a/src/MOLECULE/bond_table.h b/src/MOLECULE/bond_table.h
index dda3290d04..f0d289f055 100644
--- a/src/MOLECULE/bond_table.h
+++ b/src/MOLECULE/bond_table.h
@@ -38,7 +38,7 @@ class BondTable : public Bond {
double single(int, double, int, int);
private:
- int tabstyle,n,nm1;
+ int tabstyle,tablength;
double *r0;
struct Table {
diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp
index 1aa41e9a01..a8f2eddc50 100644
--- a/src/compute_dihedral_local.cpp
+++ b/src/compute_dihedral_local.cpp
@@ -109,9 +109,8 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
{
int i,m,n,atom1,atom2,atom3,atom4;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
- double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2;
- double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2;
- double c2mag,sin2,sc1,sc2,s1,s2,s12,c;
+ double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv;
+ double s,c;
double *pbuf;
double **x = atom->x;
@@ -148,71 +147,53 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
if (flag) {
- // phi calculation from dihedral style OPLS
+ // phi calculation from dihedral style harmonic
if (pflag >= 0) {
vb1x = x[atom1][0] - x[atom2][0];
vb1y = x[atom1][1] - x[atom2][1];
vb1z = x[atom1][2] - x[atom2][2];
domain->minimum_image(vb1x,vb1y,vb1z);
-
+
vb2x = x[atom3][0] - x[atom2][0];
vb2y = x[atom3][1] - x[atom2][1];
vb2z = x[atom3][2] - x[atom2][2];
domain->minimum_image(vb2x,vb2y,vb2z);
-
+
vb2xm = -vb2x;
vb2ym = -vb2y;
vb2zm = -vb2z;
domain->minimum_image(vb2xm,vb2ym,vb2zm);
-
+
vb3x = x[atom4][0] - x[atom3][0];
vb3y = x[atom4][1] - x[atom3][1];
vb3z = x[atom4][2] - x[atom3][2];
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;
- r12c1 = 1.0 / (b1mag*b2mag);
- c1mag = ctmp * r12c1;
+ ax = vb1y*vb2zm - vb1z*vb2ym;
+ ay = vb1z*vb2xm - vb1x*vb2zm;
+ 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;
- r12c2 = 1.0 / (b2mag*b3mag);
- c2mag = ctmp * r12c2;
-
- sin2 = MAX(1.0 - c1mag*c1mag,0.0);
- sc1 = sqrt(sin2);
- 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;
+ rasq = ax*ax + ay*ay + az*az;
+ rbsq = bx*bx + by*by + bz*bz;
+ 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);
- s1 = sc1 * sc1;
- s2 = sc2 * sc2;
- s12 = sc1 * sc2;
- c = (c0 + c1mag*c2mag) * s12;
+ c = (ax*bx + ay*by + az*bz)*rabinv;
+ s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
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;
}
diff --git a/src/pair_table.cpp b/src/pair_table.cpp
index c652a4148d..a116857108 100644
--- a/src/pair_table.cpp
+++ b/src/pair_table.cpp
@@ -76,6 +76,7 @@ void PairTable::compute(int eflag, int vflag)
Table *tb;
union_int_float_t rsq_lookup;
+ int tlm1 = tablength - 1;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@@ -127,19 +128,19 @@ void PairTable::compute(int eflag, int vflag)
if (tabstyle == LOOKUP) {
itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta);
- if (itable >= nm1)
+ if (itable >= tlm1)
error->one("Pair distance > table outer cutoff");
fpair = factor_lj * tb->f[itable];
} else if (tabstyle == LINEAR) {
itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta);
- if (itable >= nm1)
+ if (itable >= tlm1)
error->one("Pair distance > table outer cutoff");
fraction = (rsq - tb->rsq[itable]) * tb->invdelta;
value = tb->f[itable] + fraction*tb->df[itable];
fpair = factor_lj * value;
} else if (tabstyle == SPLINE) {
itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta);
- if (itable >= nm1)
+ if (itable >= tlm1)
error->one("Pair distance > table outer cutoff");
b = (rsq - tb->rsq[itable]) * tb->invdelta;
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 error->all("Unknown table style in pair_style command");
- n = force->inumeric(arg[1]);
- nm1 = n - 1;
+ tablength = force->inumeric(arg[1]);
+ if (tablength < 2) error->all("Illegal number of pair table entries");
// 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
// this is only the case if r values needed by final tables
// exactly match r values read from file
+ // for tabstyle SPLINE, always need to build spline tables
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;
- if (tabstyle == SPLINE && tb->ninput == n &&
- tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1;
- if (tabstyle == BITMAP && tb->ninput == 1 << n &&
+ if (tabstyle == BITMAP && tb->ninput == 1 << tablength &&
tb->rflag == BMP && tb->rhi == tb->cut) tb->match = 1;
-
if (tb->rflag == BMP && tb->match == 0)
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)
{
+ int tlm1 = tablength-1;
+
// inner = inner table bound
// cut = outer table bound
// 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;
else inner = tb->rfile[0];
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;
// 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
if (tabstyle == LOOKUP) {
- tb->e = (double *) memory->smalloc(nm1*sizeof(double),"pair:e");
- tb->f = (double *) memory->smalloc(nm1*sizeof(double),"pair:f");
+ tb->e = (double *) memory->smalloc(tlm1*sizeof(double),"pair:e");
+ tb->f = (double *) memory->smalloc(tlm1*sizeof(double),"pair:f");
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;
r = sqrt(rsq);
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
if (tabstyle == LINEAR) {
- tb->rsq = (double *) memory->smalloc(n*sizeof(double),"pair:rsq");
- tb->e = (double *) memory->smalloc(n*sizeof(double),"pair:e");
- tb->f = (double *) memory->smalloc(n*sizeof(double),"pair:f");
- tb->de = (double *) memory->smalloc(nm1*sizeof(double),"pair:de");
- tb->df = (double *) memory->smalloc(nm1*sizeof(double),"pair:df");
+ tb->rsq = (double *) memory->smalloc(tablength*sizeof(double),"pair:rsq");
+ tb->e = (double *) memory->smalloc(tablength*sizeof(double),"pair:e");
+ tb->f = (double *) memory->smalloc(tablength*sizeof(double),"pair:f");
+ tb->de = (double *) memory->smalloc(tlm1*sizeof(double),"pair:de");
+ tb->df = (double *) memory->smalloc(tlm1*sizeof(double),"pair:df");
double r,rsq;
- for (int i = 0; i < n; i++) {
+ for (int i = 0; i < tablength; i++) {
rsq = tb->innersq + i*tb->delta;
r = sqrt(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->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
if (tabstyle == SPLINE) {
- tb->rsq = (double *) memory->smalloc(n*sizeof(double),"pair:rsq");
- tb->e = (double *) memory->smalloc(n*sizeof(double),"pair:e");
- tb->f = (double *) memory->smalloc(n*sizeof(double),"pair:f");
- tb->e2 = (double *) memory->smalloc(n*sizeof(double),"pair:e2");
- tb->f2 = (double *) memory->smalloc(n*sizeof(double),"pair:f2");
+ tb->rsq = (double *) memory->smalloc(tablength*sizeof(double),"pair:rsq");
+ tb->e = (double *) memory->smalloc(tablength*sizeof(double),"pair:e");
+ tb->f = (double *) memory->smalloc(tablength*sizeof(double),"pair:f");
+ tb->e2 = (double *) memory->smalloc(tablength*sizeof(double),"pair:e2");
+ tb->f2 = (double *) memory->smalloc(tablength*sizeof(double),"pair:f2");
tb->deltasq6 = tb->delta*tb->delta / 6.0;
double r,rsq;
- for (int i = 0; i < n; i++) {
+ for (int i = 0; i < tablength; i++) {
rsq = tb->innersq + i*tb->delta;
r = sqrt(rsq);
tb->rsq[i] = rsq;
@@ -637,8 +638,8 @@ void PairTable::compute_table(Table *tb)
// ep0,epn = dE/dr at inner and at cut
double ep0 = - tb->f[0];
- double epn = - tb->f[nm1];
- spline(tb->rsq,tb->e,n,ep0,epn,tb->e2);
+ double epn = - tb->f[tlm1];
+ spline(tb->rsq,tb->e,tablength,ep0,epn,tb->e2);
// fp0,fpn = dh/dg at inner and at cut
// 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 =
- (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 {
double rsq2 = tb->cut * tb->cut;
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)) /
sqrt(rsq1)) / (secant_factor*tb->delta);
}
- for (int i = 0; i < n; i++) tb->f[i] /= sqrt(tb->rsq[i]);
- spline(tb->rsq,tb->f,n,fp0,fpn,tb->f2);
+ for (int i = 0; i < tablength; i++) tb->f[i] /= sqrt(tb->rsq[i]);
+ spline(tb->rsq,tb->f,tablength,fp0,fpn,tb->f2);
}
// bitmapped linear tables
@@ -683,8 +684,8 @@ void PairTable::compute_table(Table *tb)
// linear lookup tables of length ntable = 2^n
// stored value = value at lower edge of bin
- init_bitmap(inner,tb->cut,n,masklo,maskhi,tb->nmask,tb->nshiftbits);
- int ntable = 1 << n;
+ init_bitmap(inner,tb->cut,tablength,masklo,maskhi,tb->nmask,tb->nshiftbits);
+ int ntable = 1 << tablength;
int ntablem1 = ntable - 1;
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)
{
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) {
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(&n,1,MPI_INT,0,world);
- nm1 = n - 1;
+ MPI_Bcast(&tablength,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
@@ -911,23 +911,24 @@ double PairTable::single(int i, int j, int itype, int jtype, double rsq,
{
int itable;
double fraction,value,a,b,phi;
+ int tlm1 = tablength - 1;
Table *tb = &tables[tabindex[itype][jtype]];
if (rsq < tb->innersq) error->one("Pair distance < table inner cutoff");
if (tabstyle == LOOKUP) {
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];
} else if (tabstyle == LINEAR) {
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;
value = tb->f[itable] + fraction*tb->df[itable];
fforce = factor_lj * value;
} else if (tabstyle == SPLINE) {
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;
a = 1.0 - b;
value = a * tb->f[itable] + b * tb->f[itable+1] +
diff --git a/src/pair_table.h b/src/pair_table.h
index d56c1b8d40..da38a4a62d 100644
--- a/src/pair_table.h
+++ b/src/pair_table.h
@@ -40,7 +40,7 @@ class PairTable : public Pair {
void *extract(char *);
private:
- int tabstyle,n,nm1;
+ int tabstyle,tablength;
struct Table {
int ninput,rflag,fpflag,match,ntablebits;
int nshiftbits,nmask;
--
GitLab