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

Cleanup of style (removing all tabs, shortened long lines).

Browse Source
This commit is contained in:
DallasTrinkle
2017-05-04 15:28:11 -05:00
parent f7230006fe
commit 42531389df
3 changed files with 194 additions and 192 deletions
Download Patch File Download Diff File Expand all files Collapse all files

197
src/USER-MISC/pair_meam_spline.cpp
View File

@ -14,8 +14,8 @@
/* ----------------------------------------------------------------------
Contributing author: Alexander Stukowski (LLNL), alex@stukowski.com
Will Tipton (Cornell), wwt26@cornell.edu
Dallas R. Trinkle (UIUC), dtrinkle@illinois.edu
Pinchao Zhang (UIUC)
Dallas R. Trinkle (UIUC), dtrinkle@illinois.edu
Pinchao Zhang (UIUC)
see LLNL copyright notice at bottom of file
------------------------------------------------------------------------- */
@ -157,29 +157,29 @@ void PairMEAMSpline::compute(int eflag, int vflag)
double rij_sq = jdelx*jdelx + jdely*jdely + jdelz*jdelz;
if(rij_sq < cutoff*cutoff) {
double rij = sqrt(rij_sq);
double partial_sum = 0;
nextTwoBodyInfo->tag = j;
nextTwoBodyInfo->r = rij;
nextTwoBodyInfo->f = fs[i_to_potl(j)].eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->del[0] = jdelx / rij;
nextTwoBodyInfo->del[1] = jdely / rij;
nextTwoBodyInfo->del[2] = jdelz / rij;
for(int kk = 0; kk < numBonds; kk++) {
const MEAM2Body& bondk = twoBodyInfo[kk];
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
nextTwoBodyInfo->del[1]*bondk.del[1] +
nextTwoBodyInfo->del[2]*bondk.del[2]);
partial_sum += bondk.f * gs[ij_to_potl(j,bondk.tag)].eval(cos_theta);
}
rho_value += nextTwoBodyInfo->f * partial_sum;
rho_value += rhos[i_to_potl(j)].eval(rij);
numBonds++;
nextTwoBodyInfo++;
double rij = sqrt(rij_sq);
double partial_sum = 0;
nextTwoBodyInfo->tag = j;
nextTwoBodyInfo->r = rij;
nextTwoBodyInfo->f = fs[i_to_potl(j)].eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->del[0] = jdelx / rij;
nextTwoBodyInfo->del[1] = jdely / rij;
nextTwoBodyInfo->del[2] = jdelz / rij;
for(int kk = 0; kk < numBonds; kk++) {
const MEAM2Body& bondk = twoBodyInfo[kk];
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
nextTwoBodyInfo->del[1]*bondk.del[1] +
nextTwoBodyInfo->del[2]*bondk.del[2]);
partial_sum += bondk.f * gs[ij_to_potl(j,bondk.tag)].eval(cos_theta);
}
rho_value += nextTwoBodyInfo->f * partial_sum;
rho_value += rhos[i_to_potl(j)].eval(rij);
numBonds++;
nextTwoBodyInfo++;
}
}
@ -191,9 +191,9 @@ void PairMEAMSpline::compute(int eflag, int vflag)
Uprime_values[i] = Uprime_i;
if(eflag) {
if(eflag_global)
eng_vdwl += embeddingEnergy;
eng_vdwl += embeddingEnergy;
if(eflag_atom)
eatom[i] += embeddingEnergy;
eatom[i] += embeddingEnergy;
}
// Compute three-body contributions to force
@ -210,57 +210,57 @@ void PairMEAMSpline::compute(int eflag, int vflag)
MEAM2Body const* bondk = twoBodyInfo;
for(int kk = 0; kk < jj; kk++, ++bondk) {
double rik = bondk->r;
double cos_theta = (bondj.del[0]*bondk->del[0] +
bondj.del[1]*bondk->del[1] +
bondj.del[2]*bondk->del[2]);
double g_prime;
double g_value = gs[ij_to_potl(j,bondk->tag)].eval(cos_theta, g_prime);
double f_rik_prime = bondk->fprime;
double f_rik = bondk->f;
double fij = -Uprime_i * g_value * f_rik * f_rij_prime;
double fik = -Uprime_i * g_value * f_rij * f_rik_prime;
double prefactor = Uprime_i * f_rij * f_rik * g_prime;
double prefactor_ij = prefactor / rij;
double prefactor_ik = prefactor / rik;
fij += prefactor_ij * cos_theta;
fik += prefactor_ik * cos_theta;
double fj[3], fk[3];
fj[0] = bondj.del[0] * fij - bondk->del[0] * prefactor_ij;
fj[1] = bondj.del[1] * fij - bondk->del[1] * prefactor_ij;
fj[2] = bondj.del[2] * fij - bondk->del[2] * prefactor_ij;
forces_j[0] += fj[0];
forces_j[1] += fj[1];
forces_j[2] += fj[2];
fk[0] = bondk->del[0] * fik - bondj.del[0] * prefactor_ik;
fk[1] = bondk->del[1] * fik - bondj.del[1] * prefactor_ik;
fk[2] = bondk->del[2] * fik - bondj.del[2] * prefactor_ik;
forces_i[0] -= fk[0];
forces_i[1] -= fk[1];
forces_i[2] -= fk[2];
int k = bondk->tag;
forces[k][0] += fk[0];
forces[k][1] += fk[1];
forces[k][2] += fk[2];
if(evflag) {
double delta_ij[3];
double delta_ik[3];
delta_ij[0] = bondj.del[0] * rij;
delta_ij[1] = bondj.del[1] * rij;
delta_ij[2] = bondj.del[2] * rij;
delta_ik[0] = bondk->del[0] * rik;
delta_ik[1] = bondk->del[1] * rik;
delta_ik[2] = bondk->del[2] * rik;
ev_tally3(i, j, k, 0.0, 0.0, fj, fk, delta_ij, delta_ik);
}
double rik = bondk->r;
double cos_theta = (bondj.del[0]*bondk->del[0] +
bondj.del[1]*bondk->del[1] +
bondj.del[2]*bondk->del[2]);
double g_prime;
double g_value = gs[ij_to_potl(j,bondk->tag)].eval(cos_theta, g_prime);
double f_rik_prime = bondk->fprime;
double f_rik = bondk->f;
double fij = -Uprime_i * g_value * f_rik * f_rij_prime;
double fik = -Uprime_i * g_value * f_rij * f_rik_prime;
double prefactor = Uprime_i * f_rij * f_rik * g_prime;
double prefactor_ij = prefactor / rij;
double prefactor_ik = prefactor / rik;
fij += prefactor_ij * cos_theta;
fik += prefactor_ik * cos_theta;
double fj[3], fk[3];
fj[0] = bondj.del[0] * fij - bondk->del[0] * prefactor_ij;
fj[1] = bondj.del[1] * fij - bondk->del[1] * prefactor_ij;
fj[2] = bondj.del[2] * fij - bondk->del[2] * prefactor_ij;
forces_j[0] += fj[0];
forces_j[1] += fj[1];
forces_j[2] += fj[2];
fk[0] = bondk->del[0] * fik - bondj.del[0] * prefactor_ik;
fk[1] = bondk->del[1] * fik - bondj.del[1] * prefactor_ik;
fk[2] = bondk->del[2] * fik - bondj.del[2] * prefactor_ik;
forces_i[0] -= fk[0];
forces_i[1] -= fk[1];
forces_i[2] -= fk[2];
int k = bondk->tag;
forces[k][0] += fk[0];
forces[k][1] += fk[1];
forces[k][2] += fk[2];
if(evflag) {
double delta_ij[3];
double delta_ik[3];
delta_ij[0] = bondj.del[0] * rij;
delta_ij[1] = bondj.del[1] * rij;
delta_ij[2] = bondj.del[2] * rij;
delta_ik[0] = bondk->del[0] * rik;
delta_ik[1] = bondk->del[1] * rik;
delta_ik[2] = bondk->del[2] * rik;
ev_tally3(i, j, k, 0.0, 0.0, fj, fk, delta_ij, delta_ik);
}
}
forces[i][0] -= forces_j[0];
@ -304,7 +304,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
double pair_pot_deriv;
double pair_pot = phis[ij_to_potl(i,j)].eval(rij, pair_pot_deriv);
fpair += pair_pot_deriv;
fpair += pair_pot_deriv;
// Divide by r_ij to get forces from gradient
@ -406,18 +406,18 @@ void PairMEAMSpline::coeff(int narg, char **arg)
// old style: we only have one species, so we're either "NULL" or we match.
for (i = 3; i < narg; i++)
if (strcmp(arg[i],"NULL") == 0)
map[i-2] = -1;
map[i-2] = -1;
else
map[i-2] = 0;
map[i-2] = 0;
} else {
for (i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
continue;
map[i-2] = -1;
continue;
}
for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0)
break;
if (strcmp(arg[i],elements[j]) == 0)
break;
if (j < nelements) map[i-2] = j;
else error->all(FLERR,"No matching element in EAM potential file");
}
@ -460,7 +460,7 @@ void PairMEAMSpline::read_file(const char* filename)
char line[MAXLINE];
fgets(line, MAXLINE, fp);
// Second line holds potential type (currently just "meam/spline") in new potential format.
// Second line holds potential type ("meam/spline") in new potential format.
bool isNewFormat;
long loc = ftell(fp);
fgets(line, MAXLINE, fp);
@ -471,22 +471,22 @@ void PairMEAMSpline::read_file(const char* filename)
const char *sep = " ,;:-\t\n"; // overkill, but safe
word = strsep(&linep, sep);
if (! *word)
error->one(FLERR, "Need to include number of atomic species on meam/spline line in potential file");
error->one(FLERR, "Need to include number of atomic species on meam/spline line in potential file");
int n = atoi(word);
if (n<1)
error->one(FLERR, "Invalid number of atomic species on meam/spline line in potential file");
error->one(FLERR, "Invalid number of atomic species on meam/spline line in potential file");
nelements = n;
elements = new char*[n];
for (int i=0; i<n; ++i) {
word = strsep(&linep, sep);
if (! *word)
error->one(FLERR, "Not enough atomic species in meam/spline\n");
elements[i] = new char[strlen(word)+1];
strcpy(elements[i], word);
word = strsep(&linep, sep);
if (! *word)
error->one(FLERR, "Not enough atomic species in meam/spline\n");
elements[i] = new char[strlen(word)+1];
strcpy(elements[i], word);
}
} else {
isNewFormat = false;
nelements = 1; // old format only handles one species anyway; this is for backwards compatibility
nelements = 1; // old format only handles one species; (backwards compatibility)
elements = new char*[1];
elements[0] = new char[1];
strcpy(elements[0], "");
@ -607,7 +607,8 @@ double PairMEAMSpline::init_one(int i, int j)
/* ---------------------------------------------------------------------- */
int PairMEAMSpline::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
int PairMEAMSpline::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int* list_iter = list;
int* list_iter_end = list + n;
@ -646,7 +647,8 @@ double PairMEAMSpline::memory_usage()
/// Parses the spline knots from a text file.
void PairMEAMSpline::SplineFunction::parse(FILE* fp, Error* error, bool isNewFormat)
void PairMEAMSpline::SplineFunction::parse(FILE* fp, Error* error,
bool isNewFormat)
{
char line[MAXLINE];
@ -761,7 +763,8 @@ void PairMEAMSpline::SplineFunction::communicate(MPI_Comm& world, int me)
/// Writes a Gnuplot script that plots the spline function.
///
/// This function is for debugging only!
void PairMEAMSpline::SplineFunction::writeGnuplot(const char* filename, const char* title) const
void PairMEAMSpline::SplineFunction::writeGnuplot(const char* filename,
const char* title) const
{
FILE* fp = fopen(filename, "w");
fprintf(fp, "#!/usr/bin/env gnuplot\n");

160
src/USER-MISC/pair_meam_spline.h
View File

@ -28,10 +28,12 @@ PairStyle(meam/spline,PairMEAMSpline)
namespace LAMMPS_NS {
/// Set this to 1 if you intend to use MEAM potentials with non-uniform spline knots.
/// Set this to 0 if you intend to use only MEAM potentials with spline knots on a uniform grid.
///
/// With SUPPORT_NON_GRID_SPLINES == 0, the code runs about 50% faster.
// Set this to 1 if you intend to use MEAM potentials with
// non-uniform spline knots.
// Set this to 0 if you intend to use only MEAM potentials with
// spline knots on a uniform grid.
//
// With SUPPORT_NON_GRID_SPLINES == 0, the code runs about 50% faster.
#define SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES 0
@ -114,33 +116,35 @@ protected:
{
x -= xmin;
if(x <= 0.0) { // Left extrapolation.
return Y[0] + deriv0 * x;
return Y[0] + deriv0 * x;
}
else if(x >= xmax_shifted) { // Right extrapolation.
return Y[N-1] + derivN * (x - xmax_shifted);
return Y[N-1] + derivN * (x - xmax_shifted);
}
else {
#if SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES
// Do interval search.
int klo = 0;
int khi = N-1;
while(khi - klo > 1) {
int k = (khi + klo) / 2;
if(Xs[k] > x) khi = k;
else klo = k;
}
double h = Xs[khi] - Xs[klo];
// Do spline interpolation.
double a = (Xs[khi] - x)/h;
double b = 1.0 - a; // = (x - X[klo])/h
return a * Y[klo] + b * Y[khi] + ((a*a*a - a) * Y2[klo] + (b*b*b - b) * Y2[khi])*(h*h)/6.0;
// Do interval search.
int klo = 0;
int khi = N-1;
while(khi - klo > 1) {
int k = (khi + klo) / 2;
if(Xs[k] > x) khi = k;
else klo = k;
}
double h = Xs[khi] - Xs[klo];
// Do spline interpolation.
double a = (Xs[khi] - x)/h;
double b = 1.0 - a; // = (x - X[klo])/h
return a * Y[klo] + b * Y[khi] +
((a*a*a - a) * Y2[klo] + (b*b*b - b) * Y2[khi])*(h*h)/6.0;
#else
// For a spline with grid points, we can directly calculate the interval X is in.
int klo = (int)(x / h);
int khi = klo + 1;
double a = Xs[khi] - x;
double b = h - a;
return Y[khi] - a * Ydelta[klo] + ((a*a - hsq) * a * Y2[klo] + (b*b - hsq) * b * Y2[khi]);
// For a spline with regular grid, we directly calculate the interval X is in.
int klo = (int)(x / h);
int khi = klo + 1;
double a = Xs[khi] - x;
double b = h - a;
return Y[khi] - a * Ydelta[klo] +
((a*a - hsq) * a * Y2[klo] + (b*b - hsq) * b * Y2[khi]);
#endif
}
}
@ -150,37 +154,43 @@ protected:
{
x -= xmin;
if(x <= 0.0) { // Left extrapolation.
deriv = deriv0;
return Y[0] + deriv0 * x;
deriv = deriv0;
return Y[0] + deriv0 * x;
}
else if(x >= xmax_shifted) { // Right extrapolation.
deriv = derivN;
return Y[N-1] + derivN * (x - xmax_shifted);
deriv = derivN;
return Y[N-1] + derivN * (x - xmax_shifted);
}
else {
#if SPLINE_MEAM_SUPPORT_NON_GRID_SPLINES
// Do interval search.
int klo = 0;
int khi = N-1;
while(khi - klo > 1) {
int k = (khi + klo) / 2;
if(Xs[k] > x) khi = k;
else klo = k;
}
double h = Xs[khi] - Xs[klo];
// Do spline interpolation.
double a = (Xs[khi] - x)/h;
double b = 1.0 - a; // = (x - X[klo])/h
deriv = (Y[khi] - Y[klo]) / h + ((3.0*b*b - 1.0) * Y2[khi] - (3.0*a*a - 1.0) * Y2[klo]) * h / 6.0;
return a * Y[klo] + b * Y[khi] + ((a*a*a - a) * Y2[klo] + (b*b*b - b) * Y2[khi]) * (h*h) / 6.0;
// Do interval search.
int klo = 0;
int khi = N-1;
while(khi - klo > 1) {
int k = (khi + klo) / 2;
if(Xs[k] > x) khi = k;
else klo = k;
}
double h = Xs[khi] - Xs[klo];
// Do spline interpolation.
double a = (Xs[khi] - x)/h;
double b = 1.0 - a; // = (x - X[klo])/h
deriv = (Y[khi] - Y[klo]) / h +
((3.0*b*b - 1.0) * Y2[khi] -
(3.0*a*a - 1.0) * Y2[klo]) * h / 6.0;
return a * Y[klo] + b * Y[khi] +
((a*a*a - a) * Y2[klo] +
(b*b*b - b) * Y2[khi]) * (h*h) / 6.0;
#else
// For a spline with grid points, we can directly calculate the interval X is in.
int klo = (int)(x / h);
int khi = klo + 1;
double a = Xs[khi] - x;
double b = h - a;
deriv = Ydelta[klo] + ((3.0*b*b - hsq) * Y2[khi] - (3.0*a*a - hsq) * Y2[klo]);
return Y[khi] - a * Ydelta[klo] + ((a*a - hsq) * a * Y2[klo] + (b*b - hsq) * b * Y2[khi]);
// For a spline with regular grid, we directly calculate the interval X is in.
int klo = (int)(x / h);
int khi = klo + 1;
double a = Xs[khi] - x;
double b = h - a;
deriv = Ydelta[klo] + ((3.0*b*b - hsq) * Y2[khi]
- (3.0*a*a - hsq) * Y2[klo]);
return Y[khi] - a * Ydelta[klo] +
((a*a - hsq) * a * Y2[klo] + (b*b - hsq) * b * Y2[khi]);
#endif
}
}
@ -198,20 +208,20 @@ protected:
void communicate(MPI_Comm& world, int me);
private:
double* X; // Positions of spline knots
double* Xs; // Shifted positions of spline knots
double* Y; // Function values at spline knots
double* Y2; // Second derivatives at spline knots
double* Ydelta; // If this is a grid spline, Ydelta[i] = (Y[i+1]-Y[i])/h
int N; // Number of spline knots
double deriv0; // First derivative at knot 0
double derivN; // First derivative at knot (N-1)
double xmin; // The beginning of the interval on which the spline function is defined.
double xmax; // The end of the interval on which the spline function is defined.
int isGridSpline; // Indicates that all spline knots are on a regular grid.
double h; // The distance between knots if this is a grid spline with equidistant knots.
double hsq; // The squared distance between knots if this is a grid spline with equidistant knots.
double xmax_shifted; // The end of the spline interval after it has been shifted to begin at X=0.
double* X; // Positions of spline knots
double* Xs; // Shifted positions of spline knots
double* Y; // Function values at spline knots
double* Y2; // Second derivatives at spline knots
double* Ydelta; // If this is a grid spline, Ydelta[i] = (Y[i+1]-Y[i])/h
int N; // Number of spline knots
double deriv0; // First derivative at knot 0
double derivN; // First derivative at knot (N-1)
double xmin; // The beginning of the interval on which the spline function is defined.
double xmax; // The end of the interval on which the spline function is defined.
int isGridSpline;// Indicates that all spline knots are on a regular grid.
double h; // The distance between knots if this is a grid spline with equidistant knots.
double hsq; // The squared distance between knots if this is a grid spline with equidistant knots.
double xmax_shifted; // The end of the spline interval after it has been shifted to begin at X=0.
};
/// Helper data structure for potential routine.
@ -222,19 +232,19 @@ protected:
double del[3];
};
SplineFunction* phis; // Phi_i(r_ij)
SplineFunction* rhos; // Rho_ij(r_ij)
SplineFunction* fs; // f_i(r_ij)
SplineFunction* Us; // U_i(rho)
SplineFunction* gs; // g_ij(cos_theta)
double* zero_atom_energies; // Shift embedding energy by this value to make it zero for a single atom in vacuum.
SplineFunction* phis; // Phi_i(r_ij)
SplineFunction* rhos; // Rho_ij(r_ij)
SplineFunction* fs; // f_i(r_ij)
SplineFunction* Us; // U_i(rho)
SplineFunction* gs; // g_ij(cos_theta)
double* zero_atom_energies; // Shift embedding energy by this value to make it zero for a single atom in vacuum.
double cutoff; // The cutoff radius
double cutoff; // The cutoff radius
double* Uprime_values; // Used for temporary storage of U'(rho) values
int nmax; // Size of temporary array.
int maxNeighbors; // The last maximum number of neighbors a single atoms has.
MEAM2Body* twoBodyInfo; // Temporary array.
double* Uprime_values; // Used for temporary storage of U'(rho) values
int nmax; // Size of temporary array.
int maxNeighbors; // The last maximum number of neighbors a single atoms has.
MEAM2Body* twoBodyInfo; // Temporary array.
void read_file(const char* filename);
void allocate();

29
src/USER-OMP/pair_meam_spline_omp.cpp
View File

@ -140,25 +140,21 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
nextTwoBodyInfo->tag = j;
nextTwoBodyInfo->r = rij;
// nextTwoBodyInfo->f = f.eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->f = fs[i_to_potl(j)].eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->f = fs[i_to_potl(j)].eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->del[0] = jdelx / rij;
nextTwoBodyInfo->del[1] = jdely / rij;
nextTwoBodyInfo->del[2] = jdelz / rij;
for(int kk = 0; kk < numBonds; kk++) {
const MEAM2Body& bondk = myTwoBodyInfo[kk];
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
nextTwoBodyInfo->del[1]*bondk.del[1] +
nextTwoBodyInfo->del[2]*bondk.del[2]);
partial_sum += bondk.f * gs[ij_to_potl(j,bondk.tag)].eval(cos_theta);
// double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] + nextTwoBodyInfo->del[1]*bondk.del[1] + nextTwoBodyInfo->del[2]*bondk.del[2]);
// partial_sum += bondk.f * g.eval(cos_theta);
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
nextTwoBodyInfo->del[1]*bondk.del[1] +
nextTwoBodyInfo->del[2]*bondk.del[2]);
partial_sum += bondk.f * gs[ij_to_potl(j,bondk.tag)].eval(cos_theta);
}
rho_value += nextTwoBodyInfo->f * partial_sum;
// rho_value += rho.eval(rij);
rho_value += rhos[i_to_potl(j)].eval(rij);
rho_value += rhos[i_to_potl(j)].eval(rij);
numBonds++;
nextTwoBodyInfo++;
@ -167,7 +163,6 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
// Compute embedding energy and its derivative.
double Uprime_i;
// double embeddingEnergy = U.eval(rho_value, Uprime_i) - zero_atom_energy;
double embeddingEnergy = Us[i_to_potl(i)].eval(rho_value, Uprime_i)
- zero_atom_energies[i_to_potl(i)];
Uprime_thr[i] = Uprime_i;
@ -195,8 +190,7 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
+ bondj.del[1]*bondk->del[1]
+ bondj.del[2]*bondk->del[2]);
double g_prime;
// double g_value = g.eval(cos_theta, g_prime);
double g_value = gs[ij_to_potl(j,bondk->tag)].eval(cos_theta, g_prime);
double g_value = gs[ij_to_potl(j,bondk->tag)].eval(cos_theta, g_prime);
const double f_rik_prime = bondk->fprime;
const double f_rik = bondk->f;
@ -301,20 +295,15 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
if(rij_sq < cutforcesq) {
double rij = sqrt(rij_sq);
// double rho_prime;
// rho.eval(rij, rho_prime);
// double fpair = rho_prime * (Uprime_values[i] + Uprime_values[j]);
double rho_prime_i,rho_prime_j;
rhos[i_to_potl(i)].eval(rij,rho_prime_i);
rhos[i_to_potl(j)].eval(rij,rho_prime_j);
double fpair = rho_prime_j * Uprime_values[i] + rho_prime_i*Uprime_values[j];
// double pair_pot_deriv;
// double pair_pot = phi.eval(rij, pair_pot_deriv);
double pair_pot_deriv;
double pair_pot = phis[ij_to_potl(i,j)].eval(rij, pair_pot_deriv);
fpair += pair_pot_deriv;
fpair += pair_pot_deriv;
// Divide by r_ij to get forces from gradient.
fpair /= rij;