Commit 2f71245d authored by DallasTrinkle's avatar DallasTrinkle
Browse files

Removed extra "helper" functions.

parent 48072781

src/USER-MISC/pair_meam_spline.cpp

+0 −187
Original line number Diff line number Diff line
@@ -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

------------------------------------------------------------------------- */

src/USER-MISC/pair_meam_spline.h

+0 −5
Original line number Diff line number Diff line
@@ -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);

CRA Git | Maintained and supported by SUSTech CRA and CCSE