EAM Class Reference

EAM (Embedded Atom Method) potential with cell list neighbor finding. More...

#include <EAM.h>

Inheritance diagram for EAM:

Classes
struct	element_parameters

Public Member Functions
	EAM (const Parameters &params)
	~EAM () override=default
void	cleanMemory ()
void	force (long N, const double *R, const int *atomicNrs, double *F, double *U, double *variance, const double *fullbox) override
Public Member Functions inherited from eonc::Potential
	Potential (PotType a_ptype)
	Potential (PotType a_ptype, const Parameters &)
	Potential (const Parameters &a_params)
virtual	~Potential ()
std::tuple< double, AtomMatrix >	get_ef (const AtomMatrix &pos, const VectorXi &atmnrs, const Matrix3d &box)
PotType	getType () const
virtual bool	isSurrogate () const noexcept
	Whether this is a surrogate (GP) potential.
virtual bool	requiresIsolatedMoleculeLayout () const noexcept
	True for molecular QM / non-PBC backends (NWChem socket, ASE ORCA/NWChem, …).
virtual bool	isThreadSafe () const noexcept
	Whether this potential's force() can be called from multiple threads on the SAME instance.
virtual bool	isSharedInstanceThreadSafe () const noexcept
	Conservative gate for sharing one Potential instance across threads.
virtual bool	needsPerImageInstance () const noexcept
	Whether NEB should create separate Potential instances per image for true parallel force evaluation.
virtual bool	supportsBatchEvaluation () const noexcept
	Whether this potential supports batched evaluation of N systems in a single call.
virtual void	forceBatch (long nSystems, long nAtoms, const double *const *positions, const int *const *atomicNrs, double *const *forces, double *energies, double *variances, const double *const *boxes)
	Evaluate forces for N systems in a single call.

Private Member Functions
void	calc_force (long N, double *R, const int *atomicNrs, double *F, double *U, const double *box)
void	new_celllist (long N, const double *box, long *num_axis, long *cell_length, long *celllist_new, long num_cells, double *Rnew)
void	cell_to_neighbor (long N, long num_of_cells, long *num_axis, long *cell_length, long *celllist_new, long *neigh_list)
int	update_cell_list (long N, long num_cells, long *num_axis, long *cell_length, long *celllist_old, double *Rnew)
element_parameters	get_element_parameters (int atomic_number)

Static Private Member Functions
static double	embedding_function (const double *func_coeff, double rho)
static double	embedding_force (const double *func_coeff, double rho)

Private Attributes
std::vector< long >	celllist_old_
std::vector< long >	celllist_new_
std::vector< long >	neigh_list_
bool	initialized_ {false}
std::array< double, 3 >	rc_

Static Private Attributes
static const element_parameters	el_params []

Additional Inherited Members
Public Attributes inherited from eonc::Potential
std::atomic< size_t >	forceCallCounter
Protected Attributes inherited from eonc::Potential
PotType	ptype

Detailed Description

EAM (Embedded Atom Method) potential with cell list neighbor finding.

Definition at line 23 of file EAM.h.

Constructor & Destructor Documentation

EAM()

EAM::EAM ( const Parameters & params )

inlineexplicit

Definition at line 29 of file EAM.h.

      : Potential(PotType::EAM_AL, params),
        rc_{6.0, 6.0, 6.0} {}

~EAM()

EAM::~EAM ( )

overridedefault

Member Function Documentation

calc_force()

void EAM::calc_force ( long N, double * R, const int * atomicNrs, double * F, double * U, const double * box )

private

Definition at line 125 of file EAM.cpp.

                                                   {
  std::vector<double> drho_dr(3 * N);
  for (long k = 0; k < 3 * N; k++) {
    F[k] = 0.0;
  }
  *U = 0;
 
  // Pre-compute half-box for PBC
  const double halfBox0 = box[0] * 0.5;
  const double halfBox1 = box[1] * 0.5;
  const double halfBox2 = box[2] * 0.5;
 
  for (long i = 0; i < N; i++) {
    std::fill(drho_dr.begin(), drho_dr.end(), 0.0);
    element_parameters epar = get_element_parameters(atomicNrs[i]);
    const double cutoff2 = epar.r_cut * epar.r_cut;
 
    double dens = 0.0;
    const double xi = R[3 * i];
    const double yi = R[3 * i + 1];
    const double zi = R[3 * i + 2];
 
    for (long j = 0; j < N; j++) {
      if (i == j)
        continue;
 
      double dx = xi - R[3 * j];
      double dy = yi - R[3 * j + 1];
      double dz = zi - R[3 * j + 2];
 
      // Minimum image convention
      if (dx > halfBox0)
        dx -= box[0];
      else if (dx < -halfBox0)
        dx += box[0];
      if (dy > halfBox1)
        dy -= box[1];
      else if (dy < -halfBox1)
        dy += box[1];
      if (dz > halfBox2)
        dz -= box[2];
      else if (dz < -halfBox2)
        dz += box[2];
 
      double r2 = dx * dx + dy * dy + dz * dz;
      // r^2 cutoff avoids sqrt for far pairs
      if (r2 > cutoff2)
        continue;
 
      double r = std::sqrt(r2);
 
      // r^6 = (r^2)^3 instead of pow(r, 6)
      double r6 = r2 * r2 * r2;
      double exp_b1 = std::exp(-epar.beta1 * r);
      double exp_b2 = std::exp(-2.0 * epar.beta2 * r);
      double rho_pair = exp_b1 + 512.0 * exp_b2;
 
      dens += r6 * rho_pair;
 
      // r^5 = r^4 * r = (r^2)^2 * r instead of pow(r, 5)
      double r5 = r2 * r2 * r;
      double mag_force_den =
          6.0 * r5 * rho_pair +
          r6 * (-epar.beta1 * exp_b1 - 1024.0 * epar.beta2 * exp_b2);
 
      double invR = 1.0 / r;
      double fscale = -mag_force_den * invR;
      drho_dr[3 * i] += fscale * dx;
      drho_dr[3 * i + 1] += fscale * dy;
      drho_dr[3 * i + 2] += fscale * dz;
      drho_dr[3 * j] -= fscale * dx;
      drho_dr[3 * j + 1] -= fscale * dy;
      drho_dr[3 * j + 2] -= fscale * dz;
 
      if (j > i) {
        // Morse pair portion: simplify exp chain
        double expArg = std::exp(-epar.alphaM * (r - epar.Rm));
        double d = 1.0 - expArg;
        double phi_r = epar.Dm * d * d - epar.Dm;
        // Force: 2*Dm*alpha*d*(d-1) =
        // 2*Dm*alpha*exp(-a*(r-re))*(1-exp(-a*(r-re)))
        double mag_force = 2.0 * epar.alphaM * epar.Dm * d * (d - 1.0);
 
        double fcomp_scale = mag_force * invR;
        F[3 * i] -= fcomp_scale * dx;
        F[3 * i + 1] -= fcomp_scale * dy;
        F[3 * i + 2] -= fcomp_scale * dz;
        F[3 * j] += fcomp_scale * dx;
        F[3 * j + 1] += fcomp_scale * dy;
        F[3 * j + 2] += fcomp_scale * dz;
        *U += phi_r;
      }
    }
 
    *U += embedding_function(epar.func_coeff, dens);
    double dF_drho = embedding_force(epar.func_coeff, dens);
    for (long k = 0; k < 3 * N; k++) {
      F[k] -= dF_drho * drho_dr[k];
    }
  }
}

cell_to_neighbor()

void EAM::cell_to_neighbor ( long N, long num_of_cells, long * num_axis, long * cell_length, long * celllist_new, long * neigh_list )

private

Definition at line 250 of file EAM.cpp.

                                             {
  long num_cells = num_axis[0] * num_axis[1] * num_axis[2];
 
  for (long i = 0; i < N * (N + 1); i++) {
    neigh_list[i] = -1;
  }
 
  std::vector<long> neighbors(N + 1);
  std::vector<long> cell_list_copy(num_cells * (N + 1));
 
  for (long j = 0; j < num_axis[0]; j++) {
    for (long j1 = 0; j1 < num_axis[1]; j1++) {
      for (long j2 = 0; j2 < num_axis[2]; j2++) {
        std::fill(neighbors.begin(), neighbors.end(), -1);
 
        long cur_index = j * num_axis[1] * num_axis[2] + j1 * num_axis[2] + j2;
 
        std::copy(celllist_new, celllist_new + num_cells * (N + 1),
                  cell_list_copy.begin());
 
        if (cell_list_copy[cur_index * (N + 1) + N] == 0) {
          continue;
        }
 
        for (long d1 = -1; d1 < 2; d1++) {
          for (long d2 = -1; d2 < 2; d2++) {
            for (long d3 = -1; d3 < 2; d3++) {
              std::array<long, 3> pos = {j - d1, j1 - d2, j2 - d3};
              // Periodic wrapping
              for (long y = 0; y < 3; y++) {
                if (pos[y] < 0)
                  pos[y] = num_axis[y] - 1;
                else if (pos[y] >= num_axis[y])
                  pos[y] = 0;
              }
              long neigh_index = pos[0] * num_axis[1] * num_axis[2] +
                                 pos[1] * num_axis[2] + pos[2];
 
              long num = cell_list_copy[neigh_index * (N + 1) + N];
              for (long y = 0; y < num; y++) {
                long candidate = cell_list_copy[neigh_index * (N + 1) + y];
                // Check not already in neighbors
                bool already = false;
                for (long k = 0; k <= neighbors[N]; k++) {
                  if (neighbors[k] == candidate) {
                    already = true;
                    break;
                  }
                }
                if (!already) {
                  neighbors[N]++;
                  neighbors[neighbors[N]] = candidate;
                }
              }
            }
          }
        }
 
        for (long i = 0; i < celllist_new[cur_index * (N + 1) + N]; i++) {
          long cur = celllist_new[cur_index * (N + 1) + i];
          neigh_list[cur * (N + 1) + N] = 0;
          for (long temp = 0; temp <= neighbors[N]; temp++) {
            if (cur != neighbors[temp]) {
              neigh_list[cur * (N + 1) + neigh_list[cur * (N + 1) + N]] =
                  neighbors[temp];
              neigh_list[cur * (N + 1) + N]++;
            }
          }
        }
      }
    }
  }
}

cleanMemory()

void EAM::cleanMemory

(

)

Definition at line 24 of file EAM.cpp.

                      {
  // RAII: vectors clean themselves. Nothing to do.
}

embedding_force()

double EAM::embedding_force ( const double * func_coeff, double rho )

staticnodiscardprivate

Definition at line 358 of file EAM.cpp.

                                                                {
  // Derivative of embedding function via Horner's method
  double result = func_coeff[8] * 8;
  for (int i = 7; i >= 1; i--) {
    result = result * rho + i * func_coeff[i];
  }
  return -result;
}

embedding_function()

double EAM::embedding_function ( const double * func_coeff, double rho )

staticnodiscardprivate

Definition at line 349 of file EAM.cpp.

                                                                   {
  // Horner's method for 8th order polynomial
  double result = func_coeff[8];
  for (int i = 7; i >= 0; i--) {
    result = result * rho + func_coeff[i];
  }
  return result;
}

force()

void EAM::force ( long N, const double * R, const int * atomicNrs, double * F, double * U, double * variance, const double * fullbox )

overridevirtual

Implements eonc::Potential.

Definition at line 28 of file EAM.cpp.

                                                                    {
  variance = nullptr;
  std::array<double, 3> box = {fullbox[0], fullbox[4], fullbox[8]};
 
  std::array<long, 3> num_axis;
  std::array<long, 3> cell_length;
 
  for (long i = 0; i < 3; i++) {
    num_axis[i] = static_cast<long>(box[i] / rc_[i]) + 1;
  }
  for (long i = 0; i < 3; i++) {
    cell_length[i] = static_cast<long>(box[i] / (num_axis[i] - 1));
  }
 
  long num_cells = num_axis[0] * num_axis[1] * num_axis[2];
 
  if (!initialized_) {
    celllist_old_.assign(num_cells * (N + 1), 0);
    celllist_new_.assign(num_cells * (N + 1), 0);
    neigh_list_.assign(N * (N + 1), 0);
  }
 
  *U = 0;
  for (long k = 0; k < 3 * N; k++) {
    F[k] = 0.0;
  }
 
  // Find minimum coordinates for shifting all positions to positive
  double xmin = std::numeric_limits<double>::max();
  double ymin = xmin, zmin = xmin;
 
  std::vector<double> Rtemp(3 * N);
  std::vector<double> Rnew(3 * N);
 
  for (long i = 0; i < 3 * N; i += 3) {
    Rtemp[i] = R[i];
    Rtemp[i + 1] = R[i + 1];
    Rtemp[i + 2] = R[i + 2];
    xmin = std::min(xmin, R[i]);
    ymin = std::min(ymin, R[i + 1]);
    zmin = std::min(zmin, R[i + 2]);
  }
  // Only shift if coordinates are negative
  if (xmin > 0)
    xmin = 0;
  if (ymin > 0)
    ymin = 0;
  if (zmin > 0)
    zmin = 0;
 
  for (long i = 0; i < 3 * N; i += 3) {
    Rtemp[i] += std::abs(xmin);
    Rtemp[i + 1] += std::abs(ymin);
    Rtemp[i + 2] += std::abs(zmin);
  }
 
  // Enforce periodic boundary conditions
  for (long i = 0; i < 3 * N; i++) {
    while (Rtemp[i] > box[i % 3]) {
      Rtemp[i] -= box[i % 3];
    }
  }
 
  std::copy(Rtemp.begin(), Rtemp.end(), Rnew.begin());
 
  if (!initialized_) {
    new_celllist(N, box.data(), num_axis.data(), cell_length.data(),
                 celllist_new_.data(), num_cells, Rnew.data());
    cell_to_neighbor(N, num_cells, num_axis.data(), cell_length.data(),
                     celllist_new_.data(), neigh_list_.data());
  } else {
    if (update_cell_list(N, num_cells, num_axis.data(), cell_length.data(),
                         celllist_old_.data(), Rnew.data()) > 0) {
      new_celllist(N, box.data(), num_axis.data(), cell_length.data(),
                   celllist_new_.data(), num_cells, Rnew.data());
      cell_to_neighbor(N, num_cells, num_axis.data(), cell_length.data(),
                       celllist_new_.data(), neigh_list_.data());
    }
  }
 
  calc_force(N, Rnew.data(), atomicNrs, F, U, box.data());
 
  std::copy(celllist_new_.begin(), celllist_new_.end(), celllist_old_.begin());
 
  initialized_ = true;
}

get_element_parameters()

EAM::element_parameters EAM::get_element_parameters ( int atomic_number )

nodiscardprivate

Definition at line 116 of file EAM.cpp.

                                                                   {
  for (int i = 0; i < NPARAMS; i++) {
    if (el_params[i].Z == atomic_number) {
      return el_params[i];
    }
  }
  throw 14324; // Element not found
}

new_celllist()

void EAM::new_celllist ( long N, const double * box, long * num_axis, long * cell_length, long * celllist_new, long num_cells, double * Rnew )

private

Definition at line 228 of file EAM.cpp.

                                     {
  std::vector<long> cell_list(num_cells * (N + 1), -1);
  // Last slot of each cell holds atom count
  for (long i = 0; i < num_cells; i++) {
    cell_list[i * (N + 1) + N] = 0;
  }
 
  for (long i = 0; i < N; i++) {
    long cx = static_cast<long>(Rnew[3 * i] / cell_length[0]);
    long cy = static_cast<long>(Rnew[3 * i + 1] / cell_length[1]);
    long cz = static_cast<long>(Rnew[3 * i + 2] / cell_length[2]);
    long cell = cx * num_axis[1] * num_axis[2] + cy * num_axis[2] + cz;
 
    cell_list[cell * (N + 1) + cell_list[cell * (N + 1) + N]] = i;
    cell_list[cell * (N + 1) + N]++;
  }
 
  std::copy(cell_list.begin(), cell_list.end(), celllist_new);
}

update_cell_list()

int EAM::update_cell_list ( long N, long num_cells, long * num_axis, long * cell_length, long * celllist_old, double * Rnew )

private

Definition at line 326 of file EAM.cpp.

                                                                               {
  int changed = 0;
  std::vector<long> table(N);
 
  for (long i = 0; i < num_cells; i++) {
    for (long j = 0; j < celllist_old[i * (N + 1) + N]; j++) {
      table[celllist_old[i * (N + 1) + j]] = i;
    }
  }
 
  for (long i = 0; i < N; i++) {
    long cx = static_cast<long>(Rnew[3 * i] / cell_length[0]);
    long cy = static_cast<long>(Rnew[3 * i + 1] / cell_length[1]);
    long cz = static_cast<long>(Rnew[3 * i + 2] / cell_length[2]);
    long cell = cx * num_axis[1] * num_axis[2] + cy * num_axis[2] + cz;
    if (cell != table[i]) {
      changed++;
    }
  }
  return changed;
}

Member Data Documentation

celllist_new_

std::vector<long> EAM::celllist_new_

private

Definition at line 53 of file EAM.h.

celllist_old_

std::vector<long> EAM::celllist_old_

private

Definition at line 52 of file EAM.h.

el_params

const EAM::element_parameters EAM::el_params

staticprivate

Definition at line 15 of file EAM.h.

    : public Potential
#endif
{
public:
  explicit EAM(const Parameters &params)
      : Potential(PotType::EAM_AL, params),
        rc_{6.0, 6.0, 6.0} {}
 
  ~EAM() override = default;
 
  void cleanMemory();
  void force(long N, const double *R, const int *atomicNrs, double *F,
             double *U, double *variance, const double *fullbox) override;
 
private:
  struct element_parameters {
    const int Z;                // Atomic number
    const double Dm;            // Morse potential well depth
    const double alphaM;        // Curvature at Morse minimum
    const double Rm;            // Position of Morse minimum
    const double beta1;         // Density parameter 1
    const double beta2;         // Density parameter 2
    const double r_cut;         // Cutoff distance
    const double func_coeff[9]; // 8th order poly for embedding function
  };
  static const element_parameters el_params[];

initialized_

bool EAM::initialized_ {false}

private

Definition at line 55 of file EAM.h.

{false};

neigh_list_

std::vector<long> EAM::neigh_list_

private

Definition at line 54 of file EAM.h.

rc_

std::array<double, 3> EAM::rc_

private

Definition at line 56 of file EAM.h.

---

The documentation for this class was generated from the following files:

• /home/runner/work/eOn/eOn/include/eon/potentials/EAM/EAM.h
• /home/runner/work/eOn/eOn/client/potentials/EAM/EAM.cpp
• /home/runner/work/eOn/eOn/include/eon/potentials/EAM/Parameters.h