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 eonc::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
bool	isSharedInstanceThreadSafe () const noexcept override
	Cell lists live on the instance; two threads cannot share one EAM.
Public Member Functions inherited from eonc::Potential
	Potential (PotType a_ptype)
	Production default: construction-scope registry, else PotRegistry::get().
	Potential (PotType a_ptype, IPotRegistry &registry)
	Test seam: injected registry, no process-default get() counters.
	Potential (PotType a_ptype, const Parameters &p)
	Potential (const Parameters &a_params)
virtual	~Potential ()
void	force (std::span< const double > positions, std::span< const int > atomicNrs, std::span< double > forces, double energy, double variance, std::span< const double > box)
	C++ call site: size-checked view over the raw FFI force().
virtual void	setFixedMask (long nAtoms, const double *isFixed)
	Optional frozen-atom mask (nAtoms*3, 1.0 = fixed).
std::tuple< double, AtomMatrix >	get_ef (const AtomMatrix &pos, const VectorXi &atmnrs, const Matrix3d &box)
PotType	getType () const
virtual double	finiteCutoff () const noexcept
	Finite interaction range in position length units.
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 unsigned	layoutFlags () const noexcept
virtual bool	needsPerImageInstance () const noexcept
	Whether NEB should create separate Potential instances per image for true parallel force evaluation.
virtual std::shared_ptr< Potential >	clonePotential () const
	Independent instance that does not reload from disk.
virtual bool	supportsBatchEvaluation () const noexcept
	Whether this potential supports batched evaluation of N systems in a single call.
virtual bool	computesStress () const noexcept
	True when force() leaves a Cauchy stress that cauchyStress() can read until the next force() on this instance.
virtual Matrix3d	cauchyStress () const
	Cauchy stress in eV/Angstrom^3.
virtual void	forceBatchOwned (long nSystems, long nAtoms, const double const positions, const int const atomicNrs, double const forces, double energies, double variances, const double const boxes, const long *owners)
	Evaluate forces for 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)

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 Types inherited from eonc::Potential
enum class	PotLayout : unsigned { InProcess = 1u << 0 , NeedsWorkingDirectory = 1u << 1 , Subprocess = 1u << 2 }
	How the pot is executed. Combine with bitwise or. More...
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 eonc::Parameters & params )

inlineexplicit

Definition at line 29 of file EAM.h.

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

EAM::rc_

std::array< double, 3 > rc_

Definition EAM.h:60

eonc::PotType::EAM_AL

@ EAM_AL

Definition BaseStructures.h:49

◆ ~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 130 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 255 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 25 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 363 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 354 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 29 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++) {
    if (!(box[i] > 0.0) || !(rc_[i] > 0.0)) {
      throw std::invalid_argument(
          "EAM::force: box diagonal and cutoff must be positive");
    }
    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 121 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
}

◆ isSharedInstanceThreadSafe()

bool EAM::isSharedInstanceThreadSafe ( ) const

inlinenodiscardoverridevirtualnoexcept

Cell lists live on the instance; two threads cannot share one EAM.

Reimplemented from eonc::Potential.

Definition at line 39 of file EAM.h.

                                                                          {
    return false;
  }

◆ 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 233 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 331 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 57 of file EAM.h.

◆ celllist_old_

std::vector<long> EAM::celllist_old_

private

Definition at line 56 of file EAM.h.

◆ el_params

const EAM::element_parameters EAM::el_params

staticprivate

Definition at line 15 of file EAM.h.

    : public eonc::Potential
#endif
{
public:
  explicit EAM(const eonc::Parameters &params)
      : eonc::Potential(eonc::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;
  [[nodiscard]] bool isSharedInstanceThreadSafe() const noexcept override {
    return false;
  }
 
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

◆ initialized_

bool EAM::initialized_ {false}

private

Definition at line 59 of file EAM.h.

59{false};

◆ neigh_list_

std::vector<long> EAM::neigh_list_

private

Definition at line 58 of file EAM.h.

◆ rc_

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

private

Definition at line 60 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

Classes

Public Member Functions

Private Member Functions

Static Private Member Functions

Private Attributes

Static Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ EAM()

◆ ~EAM()

Member Function Documentation

◆ calc_force()

◆ cell_to_neighbor()

◆ cleanMemory()

◆ embedding_force()

◆ embedding_function()

◆ force()

◆ get_element_parameters()

◆ isSharedInstanceThreadSafe()

◆ new_celllist()

◆ update_cell_list()

Member Data Documentation

◆ celllist_new_

◆ celllist_old_

◆ el_params

◆ initialized_

◆ neigh_list_

◆ rc_