diff --git a/src/USER-MISC/pair_meam_spline.cpp b/src/USER-MISC/pair_meam_spline.cpp
index 704e76ee6118177c27b834d719f196d1ec6d9e9b..5a96d1f9c615b09169edbeccba7646cff03dd2d0 100644
--- a/src/USER-MISC/pair_meam_spline.cpp
+++ b/src/USER-MISC/pair_meam_spline.cpp
@@ -406,193 +406,6 @@ double PairMEAMSpline::three_body_density(int i)
return rho_value;
}
-double PairMEAMSpline::compute_three_body_contrib_to_charge_density(int i, int& numBonds) {
- double rho_value = 0;
- MEAM2Body* nextTwoBodyInfo = twoBodyInfo;
- double** const x = atom->x;
-
- for(int jj = 0; jj < listfull->numneigh[i]; jj++) {
- int j = listfull->firstneigh[i][jj];
- j &= NEIGHMASK;
-
- double jdelx = x[j][0] - x[i][0];
- double jdely = x[j][1] - x[i][1];
- double jdelz = x[j][2] - x[i][2];
- 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++;
- }
- }
-
- return rho_value;
-}
-
-double PairMEAMSpline::compute_embedding_energy_and_deriv(int eflag, int i, double rho_value) {
- double Uprime_i;
- double embeddingEnergy = Us[i_to_potl(i)].eval(rho_value, Uprime_i)
- - zero_atom_energies[i_to_potl(i)];
-
- Uprime_values[i] = Uprime_i;
- if(eflag) {
- if(eflag_global)
- eng_vdwl += embeddingEnergy;
- if(eflag_atom)
- eatom[i] += embeddingEnergy;
- }
- return Uprime_i;
-}
-
-void PairMEAMSpline::compute_three_body_contrib_to_forces(int i, int numBonds, double Uprime_i) {
- double forces_i[3] = {0, 0, 0};
- double** forces = atom->f;
-
- for(int jj = 0; jj < numBonds; jj++) {
- const MEAM2Body bondj = twoBodyInfo[jj];
- double rij = bondj.r;
- int j = bondj.tag;
-
- double f_rij_prime = bondj.fprime;
- double f_rij = bondj.f;
-
- double forces_j[3] = {0, 0, 0};
-
- 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);
- }
- }
-
- forces[i][0] -= forces_j[0];
- forces[i][1] -= forces_j[1];
- forces[i][2] -= forces_j[2];
- forces[j][0] += forces_j[0];
- forces[j][1] += forces_j[1];
- forces[j][2] += forces_j[2];
- }
-
- forces[i][0] += forces_i[0];
- forces[i][1] += forces_i[1];
- forces[i][2] += forces_i[2];
-}
-
-void PairMEAMSpline::compute_two_body_pair_interactions() {
- double** const x = atom->x;
- double** forces = atom->f;
-
- for(int ii = 0; ii < listhalf->inum; ii++) {
- int i = listhalf->ilist[ii];
-
- for(int jj = 0; jj < listhalf->numneigh[i]; jj++) {
- int j = listhalf->firstneigh[i][jj];
- j &= NEIGHMASK;
-
- double jdel[3];
- jdel[0] = x[j][0] - x[i][0];
- jdel[1] = x[j][1] - x[i][1];
- jdel[2] = x[j][2] - x[i][2];
- double rij_sq = jdel[0]*jdel[0] + jdel[1]*jdel[1] + jdel[2]*jdel[2];
-
- if(rij_sq < cutoff*cutoff) {
- double rij = sqrt(rij_sq);
-
- 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 = phis[ij_to_potl(i,j)].eval(rij, pair_pot_deriv);
-
- fpair += pair_pot_deriv;
-
- // Divide by r_ij to get forces from gradient
-
- fpair /= rij;
-
- forces[i][0] += jdel[0]*fpair;
- forces[i][1] += jdel[1]*fpair;
- forces[i][2] += jdel[2]*fpair;
- forces[j][0] -= jdel[0]*fpair;
- forces[j][1] -= jdel[1]*fpair;
- forces[j][2] -= jdel[2]*fpair;
- if (evflag) ev_tally(i, j, atom->nlocal, force->newton_pair,
- pair_pot, 0.0, -fpair, jdel[0], jdel[1], jdel[2]);
- }
- }
- }
-}
-
/* ----------------------------------------------------------------------
helper functions to map atom types to potential array indices
------------------------------------------------------------------------- */
diff --git a/src/USER-MISC/pair_meam_spline.h b/src/USER-MISC/pair_meam_spline.h
index 7e4f421909be04e3fe6f1cb837e56e06b4f99234..4be8f81ede732fc8d6146737a0c3c4b8457300c0 100644
--- a/src/USER-MISC/pair_meam_spline.h
+++ b/src/USER-MISC/pair_meam_spline.h
@@ -52,11 +52,6 @@ public:
// helper functions for compute()
- double compute_three_body_contrib_to_charge_density(int i, int& numBonds); // returns rho_value and returns numBonds by reference
- double compute_embedding_energy_and_deriv(int eflag, int i, double rho_value); // returns the derivative of the embedding energy Uprime_i
- void compute_three_body_contrib_to_forces(int i, int numBonds, double Uprime_i);
- void compute_two_body_pair_interactions();
-
int ij_to_potl(int i, int j);
int i_to_potl(int i);