LAMMPSPot Class Reference

#include <LAMMPSPot.h>

Inheritance diagram for LAMMPSPot:

Public Member Functions
bool	needsPerImageInstance () const noexcept override
	Whether NEB should create separate Potential instances per image for true parallel force evaluation.
	LAMMPSPot (const eonc::Parameters &p)
	Production: process-default LammpsLoader and POSIX worker isolation.
	LAMMPSPot (const eonc::Parameters &p, eonc::ILammpsLoader &loader)
	Test seam: injected loader, no worker fork, no process-default load.
	~LAMMPSPot ()
void	cleanMemory ()
void	force (long N, const double *R, const int *atomicNrs, double *F, double *U, double *variance, const double *box) override
void	setFixedMask (long nAtoms, const double *isFixed) override
	Optional frozen-atom mask (nAtoms*3, 1.0 = fixed).
bool	computesStress () const noexcept override
	True when force() leaves a Cauchy stress that cauchyStress() can read until the next force() on this instance.
Matrix3d	cauchyStress () const override
	Cauchy stress in eV/Angstrom^3.
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().
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	isSharedInstanceThreadSafe () const noexcept
	Conservative gate for sharing one Potential instance across threads.
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 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
	LAMMPSPot (const eonc::Parameters &p, eonc::ILammpsLoader &loader, bool isolate_worker)
void	makeNewLAMMPS (long N, const double *R, const int *atomicNrs, const double *box)
void	applySetforce (long N)
void	ensureWorker ()
bool	noteScreenGeometry (long N, const double *box)
void	runWorkerLoop ()
void	stopWorker ()
void	forceLocal (long N, const double *R, const int *atomicNrs, double *F, double *U, const double *box)
void	lammpsCommand (const char *cmd)
void	drainLammpsScreen ()

Private Attributes
eonc::ILammpsLoader &	loader_
int	lammpsThr {0}
bool	lammpsLogging_ {false}
int	lammpsLogIndex_ {0}
std::mutex	maskMutex_
std::string	lammpsScreenPath_
std::int64_t	lammpsScreenPos_ {0}
bool	lammpsScreenRestart_ {false}
bool	workerChild_ {false}
long	numberOfAtoms {0}
double	oldBox [9] {}
void *	LAMMPSObj {nullptr}
bool	realunits {false}
Matrix3d	stress_ {Matrix3d::Zero()}
bool	haveStress_ {false}
std::vector< double >	fixedMask_
long	maskN_ {0}
std::mutex	workerMutex
int	workerRespawnsLeft {3}
int	workerPid {-1}
int	reqFd {-1}
int	resFd {-1}
bool	workerSpawned {false}
bool	screenHaveGeom_ {false}
long	screenAtoms_ {0}
double	screenBox_ [9] {}

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

Definition at line 70 of file LAMMPSPot.h.

Constructor & Destructor Documentation

LAMMPSPot() [1/3]

LAMMPSPot::LAMMPSPot ( const eonc::Parameters & p )

explicit

Production: process-default LammpsLoader and POSIX worker isolation.

Definition at line 44 of file LAMMPSPot.cpp.

    : LAMMPSPot(p, eonc::LammpsLoader::instance(), true) {}

LAMMPSPot() [2/3]

LAMMPSPot::LAMMPSPot	(	const eonc::Parameters &	p,
		eonc::ILammpsLoader &	loader )

Test seam: injected loader, no worker fork, no process-default load.

Definition at line 47 of file LAMMPSPot.cpp.

    : LAMMPSPot(p, loader, false) {}

~LAMMPSPot()

LAMMPSPot::~LAMMPSPot

(

)

Definition at line 85 of file LAMMPSPot.cpp.

{ cleanMemory(); }

LAMMPSPot() [3/3]

LAMMPSPot::LAMMPSPot ( const eonc::Parameters & p, eonc::ILammpsLoader & loader, bool isolate_worker )

private

Definition at line 50 of file LAMMPSPot.cpp.

    : eonc::Potential(p),
      loader_{loader},
      lammpsThr{p.potential_options().LAMMPSThreads},
      lammpsLogging_{p.potential_options().LAMMPSLogging}
#ifdef EONMPI
      ,
      mpiComm{eonc::getMpiClientComm(p)}
#endif
{
  // Fail fast if LAMMPS library not available
  loader_.require_loaded();
  if (lammpsLogging_) {
    static std::atomic<int> ids{0};
    lammpsLogIndex_ = ids.fetch_add(1);
    lammpsScreenPath_ =
        "client_lammps-" + std::to_string(lammpsLogIndex_) + ".log";
  }
#if !defined(EONMPI) && !defined(IS_WINDOWS)
  if (!isolate_worker) {
    return;
  }
  // Fork the worker NOW, at construction, before this process ever opens a
  // LAMMPS instance (and thus before liblammps initialises MPI).  Open MPI
  // does not support using MPI in a process that called MPI_Init before fork,
  // so the worker must be spawned from a still-MPI-clean parent.  Every
  // LAMMPSPot -- endpoints and per-image alike -- runs its LAMMPS in its own
  // child process, so the parent never initialises MPI at all.
  ensureWorker();
#else
  (void)isolate_worker;
#endif
}

Member Function Documentation

applySetforce()

void LAMMPSPot::applySetforce ( long N )

private

Definition at line 101 of file LAMMPSPot.cpp.

                                    {
  std::vector<double> mask;
  {
    // The worker child inherits this mutex across fork. Only the parent locks.
    std::unique_lock<std::mutex> lock(maskMutex_, std::defer_lock);
    if (!workerChild_) {
      lock.lock();
    }
    if (LAMMPSObj == nullptr || maskN_ != N || fixedMask_.empty()) {
      return;
    }
    mask = fixedMask_;
  }
  static constexpr const char *kUnfix[] = {"unfix eon_fx", "unfix eon_fy",
                                           "unfix eon_fz", "unfix eon_freeze"};
  static constexpr const char *kUngroup[] = {
      "group eon_fx delete", "group eon_fy delete", "group eon_fz delete",
      "group eon_frozen delete"};
  for (const char *cmd : kUnfix) {
    try {
      lammpsCommand(cmd);
    } catch (...) {
    }
  }
  for (const char *cmd : kUngroup) {
    try {
      lammpsCommand(cmd);
    } catch (...) {
    }
  }
  // Matter stores per-axis isFixed. A z-only freeze must not leave LAMMPS
  // free to move that atom in z.
  std::string ids[3];
  for (long i = 0; i < N; ++i) {
    for (int ax = 0; ax < 3; ++ax) {
      if (mask[static_cast<size_t>(3 * i + ax)] >= 0.5) {
        ids[ax] += std::format("{} ", i + 1);
      }
    }
  }
  static constexpr const char *kGroup[] = {"eon_fx", "eon_fy", "eon_fz"};
  static constexpr const char *kFix[] = {
      "fix eon_fx eon_fx setforce 0.0 NULL NULL",
      "fix eon_fy eon_fy setforce NULL 0.0 NULL",
      "fix eon_fz eon_fz setforce NULL NULL 0.0"};
  for (int ax = 0; ax < 3; ++ax) {
    if (ids[ax].empty()) {
      continue;
    }
    lammpsCommand(
        ("group " + std::string(kGroup[ax]) + " id " + ids[ax]).c_str());
    lammpsCommand(kFix[ax]);
  }
}

cauchyStress()

Matrix3d LAMMPSPot::cauchyStress ( ) const

inlinenodiscardoverridevirtual

Cauchy stress in eV/Angstrom^3.

sigma = (1/V) dE/dε for the right strain h <- h (I+ε) at fixed fractional coordinates.

Reimplemented from eonc::Potential.

Definition at line 90 of file LAMMPSPot.h.

{ return stress_; }

cleanMemory()

void LAMMPSPot::cleanMemory

(

)

Definition at line 156 of file LAMMPSPot.cpp.

                            {
#if !defined(EONMPI) && !defined(IS_WINDOWS)
  stopWorker();
#endif
  if (LAMMPSObj != nullptr) {
    drainLammpsScreen();
    loader_.close(LAMMPSObj);
    LAMMPSObj = nullptr;
  }
}

computesStress()

bool LAMMPSPot::computesStress ( ) const

inlinenodiscardoverridevirtualnoexcept

True when force() leaves a Cauchy stress that cauchyStress() can read until the next force() on this instance.

Reimplemented from eonc::Potential.

Definition at line 86 of file LAMMPSPot.h.

                                                              {
    return haveStress_;
  }

drainLammpsScreen()

void LAMMPSPot::drainLammpsScreen ( )

private

Definition at line 687 of file LAMMPSPot.cpp.

                                  {
  // makeNewLAMMPS runs in the worker. The child writes the screen file and
  // must not call the process logger; the parent copies it after the reply.
  if (workerChild_ || !lammpsLogging_ || lammpsScreenPath_.empty()) {
    return;
  }
  std::error_code ec;
  const auto sz = std::filesystem::file_size(lammpsScreenPath_, ec);
  if (ec) {
    return;
  }
  const auto fileSize = static_cast<std::int64_t>(sz);
  lammpsScreenPos_ = eonc::lammpsScreenCursor(lammpsScreenPos_, fileSize,
                                              lammpsScreenRestart_);
  lammpsScreenRestart_ = false;
  std::ifstream in(lammpsScreenPath_, std::ios::in | std::ios::binary);
  if (!in) {
    return;
  }
  in.seekg(static_cast<std::streamoff>(lammpsScreenPos_));
  std::string line;
  while (std::getline(in, line)) {
    if (!line.empty() && line.back() == '\r') {
      line.pop_back();
    }
    if (!line.empty()) {
      EONC_LOG_INFO("{}", line);
    }
  }
  in.clear();
  in.seekg(0, std::ios::end);
  const auto end = in.tellg();
  if (end >= 0) {
    lammpsScreenPos_ = static_cast<std::int64_t>(end);
  }
}

ensureWorker()

void LAMMPSPot::ensureWorker ( )

private

Definition at line 234 of file LAMMPSPot.cpp.

                             {
  if (workerSpawned)
    return;
 
  int reqPipe[2]; // parent -> child
  int resPipe[2]; // child -> parent
  if (pipe(reqPipe) != 0) {
    throw std::runtime_error("LAMMPSPot: failed to create worker pipes");
  }
  if (pipe(resPipe) != 0) {
    close(reqPipe[0]);
    close(reqPipe[1]);
    throw std::runtime_error("LAMMPSPot: failed to create worker pipes");
  }
 
  // Fork BEFORE opening any LAMMPS instance in this process, so MPI is first
  // initialised inside the child.  Each child is its own process with its own
  // MPI_COMM_WORLD; concurrent children never share a communicator.
  pid_t pid = fork();
  if (pid < 0) {
    close(reqPipe[0]);
    close(reqPipe[1]);
    close(resPipe[0]);
    close(resPipe[1]);
    throw std::runtime_error("LAMMPSPot: fork for worker failed");
  }
 
  if (pid == 0) {
    // Child: keep reqPipe read end and resPipe write end.
    close(reqPipe[1]);
    close(resPipe[0]);
    reqFd = reqPipe[0];
    resFd = resPipe[1];
    // Client main arms feenableexcept unconditionally. The worker inherits
    // that mask; LAMMPS EAM (PairEAM::compute) performs IEEE divisions that
    // can raise FE_DIVBYZERO on near-coincident pairs during saddle / product
    // minimisations. Under trapping that SIGFPEs, and the continue handler
    // without MXCSR masking re-storms forever in the child (GB of identical
    // "FPE (continuing)" lines). External pot code expects soft IEEE defaults;
    // demote trapping for the whole worker process before any LAMMPS call.
    eonc::disableFPE();
    runWorkerLoop(); // never returns
  }
 
  // Parent: keep reqPipe write end and resPipe read end.
  //
  // Writing to a worker that has already exited raises SIGPIPE, whose default
  // action kills the client outright -- before writeExact can return the
  // error the caller is written to handle. A search that had converged on a
  // saddle at 0.062 eV died this way with signal 13 as its endpoints were
  // about to be minimised. Ignoring it turns the same condition into an
  // EPIPE return, which reaches the geometry-rejection path and respawns.
  std::signal(SIGPIPE, SIG_IGN);
  close(reqPipe[0]);
  close(resPipe[1]);
  reqFd = reqPipe[1];
  resFd = resPipe[0];
  workerPid = pid;
  workerSpawned = true;
}

force()

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

overridevirtual

Implements eonc::Potential.

Definition at line 412 of file LAMMPSPot.cpp.

                                                                      {
  variance = nullptr;
 
#ifdef EONMPI
  std::lock_guard<std::mutex> stateLock(workerMutex);
  forceLocal(N, R, atomicNrs, F, U, box);
#elif defined(IS_WINDOWS)
  // No fork/pipe on Windows; call forceLocal directly.
  std::lock_guard<std::mutex> stateLock(workerMutex);
  forceLocal(N, R, atomicNrs, F, U, box);
#else
  // Drive the dedicated worker process so this image's LAMMPS runs in its own
  // process (own MPI_COMM_WORLD).  Per-image NEB threads thus evaluate forces
  // as truly concurrent processes with no shared-communicator contention.
  std::lock_guard<std::mutex> workerLock(workerMutex);
  if (workerRespawnsLeft <= 0) {
    rejectGeometry(U, F, N);
    return;
  }
  const bool newWorker = !workerSpawned;
  ensureWorker();
  const bool screenRestart = newWorker || noteScreenGeometry(N, box);
 
  std::vector<double> mask(static_cast<size_t>(3 * N), 0.0);
  {
    std::lock_guard<std::mutex> lock(maskMutex_);
    if (maskN_ == N && fixedMask_.size() == static_cast<size_t>(3 * N)) {
      mask = fixedMask_;
    }
  }
  if (!writeExact(reqFd, &N, sizeof(N)) ||
      !writeExact(reqFd, atomicNrs, sizeof(int) * static_cast<size_t>(N)) ||
      !writeExact(reqFd, box, sizeof(double) * 9) ||
      !writeExact(reqFd, R, sizeof(double) * static_cast<size_t>(3 * N)) ||
      !writeExact(reqFd, mask.data(),
                  sizeof(double) * static_cast<size_t>(3 * N))) {
    // A worker stopped by an earlier rejected geometry leaves the request
    // pipe closed, so the first send after it fails. Respawning happens on
    // the next evaluation; reject this one rather than end the client.
    --workerRespawnsLeft;
    EONC_LOG_WARNING("[LAMMPSPot] send to worker failed; {} respawns left "
                     "(eon-7416)",
                     workerRespawnsLeft);
    stopWorker();
    rejectGeometry(U, F, N);
    return;
  }
 
  // eon-7416: a worker stuck on a pathological geometry (LAMMPS spinning, or
  // a NaN it never returns) would make the blocking read below hang until the
  // akmc pass times out. Bound the wait: if the worker is silent past a
  // generous per-eval deadline, kill and reap it (a stuck worker never reaches
  // EOF, so a plain waitpid would block too) and fail this evaluation so the
  // search discards the geometry; ensureWorker respawns on the next call.
  {
    struct pollfd pfd;
    pfd.fd = resFd;
    pfd.events = POLLIN;
    pfd.revents = 0;
    int pr = poll(&pfd, 1, 90000); // 90 s: orders beyond a normal force eval
    if (pr <= 0) {
      if (workerPid > 0) {
        kill(workerPid, SIGKILL);
      }
      --workerRespawnsLeft;
      EONC_LOG_WARNING(
          "[LAMMPSPot] worker force eval timed out; {} respawns left "
          "(eon-7416)",
          workerRespawnsLeft);
      stopWorker();
      rejectGeometry(U, F, N);
      return;
    }
  }
  int status = 0;
  double stressRaw[9] = {};
  if (!readExact(resFd, &status, sizeof(status)) ||
      !readExact(resFd, U, sizeof(double)) ||
      !readExact(resFd, F, sizeof(double) * static_cast<size_t>(3 * N)) ||
      !readExact(resFd, stressRaw, sizeof(stressRaw))) {
    --workerRespawnsLeft;
    EONC_LOG_WARNING(
        "[LAMMPSPot] worker died during force eval; {} respawns left "
        "(eon-7416)",
        workerRespawnsLeft);
    stopWorker();
    rejectGeometry(U, F, N);
    return;
  }
  if (screenRestart) {
    lammpsScreenRestart_ = true;
  }
  stress_.setZero();
  for (int row = 0; row < 3; ++row) {
    for (int col = 0; col < 3; ++col) {
      stress_(row, col) = stressRaw[row * 3 + col];
    }
  }
  haveStress_ = status == 0 && stress_.allFinite();
  if (status != 0) {
    drainLammpsScreen();
    --workerRespawnsLeft;
    EONC_LOG_WARNING(
        "[LAMMPSPot] worker reported an evaluation error; {} respawns left "
        "(eon-7416)",
        workerRespawnsLeft);
    stopWorker();
    rejectGeometry(U, F, N);
    return;
  }
  drainLammpsScreen();
  // A saddle search that never terminates silently truncates the event
  // table: the KMC residence time is 1/sum_j k_j over the discovered
  // mechanisms, so a dropped search removes a term and biases the clock
  // (Pedersen and Jónsson, Math. Comput. Simul. 80, 1487 (2010),
  // doi:10.1016/j.matcom.2009.02.010, Fig. 1; Alexander and Schuh,
  // Modelling Simul. Mater. Sci. Eng. 24, 065014 (2016),
  // doi:10.1088/0965-0393/24/6/065014, on catalog completeness).
  // An over-aggressive saddle-search kick can drive atoms on top of
  // each other; the EAM force overflows to NaN/Inf and LAMMPS returns it
  // rather than crashing. The min-mode search then spins on non-finite
  // gradients until the akmc pass times out (0 processes). Reject the
  // evaluation so the search discards that displacement and continues.
  // Reject the geometry, not the process. Nothing between this call and
  // main catches a throw here, so the client terminates: a search that was
  // making progress is lost, and every other search sharing the pass goes
  // with it. A large finite energy with zeroed forces reads to the
  // optimiser as an impassable wall, so it backs out of the step and the
  // search abandons this configuration and carries on.
  bool nonfinite = !std::isfinite(*U);
  for (long i = 0; i < 3 * N && !nonfinite; ++i) {
    nonfinite = !std::isfinite(F[i]);
  }
  if (nonfinite) {
    EONC_LOG_WARNING("[LAMMPSPot] non-finite force or energy; rejecting "
                     "geometry (eon-7416)");
    rejectGeometry(U, F, N);
  }
#endif
}

forceLocal()

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

private

Definition at line 554 of file LAMMPSPot.cpp.

                                                                    {
  // Same contract as ASE / Metatomic: external pot libraries are not written
  // for FE traps. Cover in-process (EONMPI / Windows) and any path that still
  // has trapping armed when forceLocal runs. Always restore so FE traps do not
  // stay demoted for the rest of the process after the first force call.
  eonc::FPEHandler fpeh;
  fpeh.eat_fpe();
  try {
    auto &lmp = loader_;
 
    bool newLammps = false;
    for (int i = 0; i < 9; i++) {
      if (oldBox[i] != box[i])
        newLammps = true;
    }
    if (numberOfAtoms != N)
      newLammps = true;
    if (newLammps) {
      makeNewLAMMPS(N, R, atomicNrs, box);
    }
    if (!LAMMPSObj) {
      throw std::runtime_error("Should have a LAMMPS instance by now");
    }
 
    lmp.scatter_atoms(LAMMPSObj, "x", 1, 3, const_cast<double *>(R));
    applySetforce(N);
    // New instance / box change: rebuild neighbors. create_atoms sits at
    // the origin; pre no would evaluate on that neighbor list.
    if (newLammps) {
      lammpsCommand("run 1 pre yes post no");
    } else {
      lammpsCommand("run 1 pre no post no");
    }
 
    auto *pe =
        static_cast<double *>(lmp.extract_variable(LAMMPSObj, "pe", nullptr));
    auto *fx =
        static_cast<double *>(lmp.extract_variable(LAMMPSObj, "fx", "all"));
    auto *fy =
        static_cast<double *>(lmp.extract_variable(LAMMPSObj, "fy", "all"));
    auto *fz =
        static_cast<double *>(lmp.extract_variable(LAMMPSObj, "fz", "all"));
    if (!pe || !fx || !fy || !fz) {
      if (pe)
        free(pe);
      if (fx)
        free(fx);
      if (fy)
        free(fy);
      if (fz)
        free(fz);
      throw std::runtime_error(
          "LAMMPS: extract_variable returned null (pe/fx/fy/fz)");
    }
    *U = *pe;
    free(pe);
 
    for (long i = 0; i < N; i++) {
      F[3 * i + 0] = fx[i];
      F[3 * i + 1] = fy[i];
      F[3 * i + 2] = fz[i];
    }
 
    // Convert kCal/mol -> eV if LAMMPS is using real units
    if (realunits) {
      constexpr double kcalPerEv = 23.0609;
      *U /= kcalPerEv;
      for (long i = 0; i < 3 * N; i++) {
        F[i] /= kcalPerEv;
      }
    }
 
    free(fx);
    free(fy);
    free(fz);
 
    // LAMMPS pressure is positive when the cell pushes outward. sigma =
    // (1/V) dE/dε has the opposite sign. metal reports bars, real reports
    // atmospheres. 1 eV/Angstrom^3 = 1.6021766208e6 bar.
    auto *sxx = static_cast<double *>(
        lmp.extract_variable(LAMMPSObj, "eon_sxx", nullptr));
    auto *syy = static_cast<double *>(
        lmp.extract_variable(LAMMPSObj, "eon_syy", nullptr));
    auto *szz = static_cast<double *>(
        lmp.extract_variable(LAMMPSObj, "eon_szz", nullptr));
    auto *sxy = static_cast<double *>(
        lmp.extract_variable(LAMMPSObj, "eon_sxy", nullptr));
    auto *sxz = static_cast<double *>(
        lmp.extract_variable(LAMMPSObj, "eon_sxz", nullptr));
    auto *syz = static_cast<double *>(
        lmp.extract_variable(LAMMPSObj, "eon_syz", nullptr));
    const bool got = sxx && syy && szz && sxy && sxz && syz;
    auto release = [&]() {
      if (sxx)
        free(sxx);
      if (syy)
        free(syy);
      if (szz)
        free(szz);
      if (sxy)
        free(sxy);
      if (sxz)
        free(sxz);
      if (syz)
        free(syz);
    };
    if (!got) {
      release();
      haveStress_ = false;
      stress_.setZero();
    } else {
      constexpr double barPerEvA3 = 1.6021766208e6;
      constexpr double barPerAtm = 1.01325;
      const double toEv =
          realunits ? -(barPerAtm / barPerEvA3) : -(1.0 / barPerEvA3);
      stress_.setZero();
      stress_(0, 0) = (*sxx) * toEv;
      stress_(1, 1) = (*syy) * toEv;
      stress_(2, 2) = (*szz) * toEv;
      stress_(0, 1) = stress_(1, 0) = (*sxy) * toEv;
      stress_(0, 2) = stress_(2, 0) = (*sxz) * toEv;
      stress_(1, 2) = stress_(2, 1) = (*syz) * toEv;
      haveStress_ = stress_.allFinite();
      release();
    }
  } catch (...) {
    fpeh.restore_fpe();
    throw;
  }
  fpeh.restore_fpe();
}

lammpsCommand()

void LAMMPSPot::lammpsCommand ( const char * cmd )

private

Definition at line 724 of file LAMMPSPot.cpp.

                                             {
  loader_.command(LAMMPSObj, cmd);
  // The forked worker must not call the process logger. The parent copies
  // the screen file after the child returns the force.
  if (!workerChild_) {
    drainLammpsScreen();
  }
}

makeNewLAMMPS()

void LAMMPSPot::makeNewLAMMPS ( long N, const double * R, const int * atomicNrs, const double * box )

private

Definition at line 733 of file LAMMPSPot.cpp.

                                                 {
  auto &lmp = loader_;
 
  numberOfAtoms = N;
  std::memcpy(oldBox, box, 9 * sizeof(double));
 
  if (LAMMPSObj != nullptr) {
    drainLammpsScreen();
    loader_.close(LAMMPSObj);
    LAMMPSObj = nullptr;
  }
 
  // Map atomic numbers to LAMMPS type indices (1-based)
  std::map<int, int> type_map;
  int ntypes = 0;
  for (long i = 0; i < N; i++) {
    if (type_map.count(atomicNrs[i]) == 0) {
      type_map.insert({atomicNrs[i], ++ntypes});
    }
  }
 
#ifdef EONMPI
  const std::vector<std::string> argStore =
      eonc::lammpsOpenArgs(lammpsLogging_, true, lammpsScreenPath_);
  std::vector<char *> argPtrs;
  argPtrs.reserve(argStore.size());
  for (const std::string &arg : argStore) {
    argPtrs.push_back(const_cast<char *>(arg.c_str()));
  }
  int lmpargc = static_cast<int>(argPtrs.size());
  char **lmpargv = argPtrs.data();
  if (!lmp.open_mpi) {
    throw std::runtime_error(
        "LAMMPS library found but lacks MPI support (lammps_open not found).\n"
        "Install an MPI-enabled LAMMPS build.");
  }
  MPI_Comm inst_comm = MPI_COMM_NULL;
  MPI_Comm_dup(mpiComm, &inst_comm); // private comm per per-image instance
  LAMMPSObj = lmp.open_mpi(lmpargc, lmpargv, inst_comm, nullptr);
#else
  const std::vector<std::string> argStore =
      eonc::lammpsOpenArgs(lammpsLogging_, false, lammpsScreenPath_);
  std::vector<char *> argPtrs;
  argPtrs.reserve(argStore.size());
  for (const std::string &arg : argStore) {
    argPtrs.push_back(const_cast<char *>(arg.c_str()));
  }
  int lmpargc = static_cast<int>(argPtrs.size());
  LAMMPSObj = lmp.open_no_mpi(lmpargc, argPtrs.data(), nullptr);
#endif
  // -screen opens with truncation. The next copy starts at the new file.
  lammpsScreenRestart_ = true;
 
  if (lammpsThr > 0) {
    std::string cmd = std::format("package omp {} force/neigh", lammpsThr);
    lammpsCommand(cmd.c_str());
  }
 
  // Detect units from in.lammps: look for "#!units real" marker
  realunits = false;
  if (std::filesystem::exists("in.lammps")) {
    std::ifstream infile("in.lammps");
    std::string line;
    while (std::getline(infile, line)) {
      if (line == "#!units real") {
        realunits = true;
        break;
      }
    }
  } else {
    if (LAMMPSObj != nullptr) {
      drainLammpsScreen();
      lmp.close(LAMMPSObj);
      LAMMPSObj = nullptr;
    }
    throw std::runtime_error(
        "LAMMPS: in.lammps not found in working directory");
  }
 
  if (realunits) {
    lammpsCommand("units real");
  } else {
    lammpsCommand("units metal");
  }
 
  lammpsCommand("atom_style charge");
  lammpsCommand("atom_modify map array sort 0 0");
  lammpsCommand("neigh_modify delay 1");
 
  // LAMMPS restricted triclinic: (ax, by, cz, bx, cx, cy).
  // Row-major Matter cell also has ay, az, bz at box[1], box[2], box[5].
  constexpr double kTiltCut = 1.0e-8;
  if (std::abs(box[1]) > kTiltCut || std::abs(box[2]) > kTiltCut ||
      std::abs(box[5]) > kTiltCut) {
    throw std::runtime_error(
        "LAMMPS: cell is not restricted triclinic (ay/az/bz must be ~0)");
  }
  std::string region_cmd =
      std::format("region cell prism 0 {} 0 {} 0 {} {} {} {} units box", box[0],
                  box[4], box[8], box[3], box[6], box[7]);
  lammpsCommand(region_cmd.c_str());
 
  std::string create_box_cmd = std::format("create_box {} cell", ntypes);
  lammpsCommand(create_box_cmd.c_str());
 
  // Initialize atoms
  for (long i = 0; i < N; i++) {
    std::string atom_cmd =
        std::format("create_atoms {} single {} {} {} units box",
                    type_map[atomicNrs[i]], 0.0, 0.0, 0.0);
    lammpsCommand(atom_cmd.c_str());
  }
 
  lammpsCommand("mass * 1.0");
 
  // Load user LAMMPS input script
  lmp.file(LAMMPSObj, "in.lammps");
 
  // Define variables for force/energy extraction
  lammpsCommand("variable fx atom fx");
  lammpsCommand("variable fy atom fy");
  lammpsCommand("variable fz atom fz");
  lammpsCommand("variable pe equal pe");
  // Kinetic part omitted. The six components are pxx pyy pzz pxy pxz pyz.
  lammpsCommand("compute eon_press all pressure NULL virial");
  lammpsCommand("variable eon_sxx equal c_eon_press[1]");
  lammpsCommand("variable eon_syy equal c_eon_press[2]");
  lammpsCommand("variable eon_szz equal c_eon_press[3]");
  lammpsCommand("variable eon_sxy equal c_eon_press[4]");
  lammpsCommand("variable eon_sxz equal c_eon_press[5]");
  lammpsCommand("variable eon_syz equal c_eon_press[6]");
}

needsPerImageInstance()

bool LAMMPSPot::needsPerImageInstance ( ) const

inlinenodiscardoverridevirtualnoexcept

Whether NEB should create separate Potential instances per image for true parallel force evaluation.

When true, NEB calls makePotential() once per image instead of sharing one instance. Override in potentials that use internal mutexes (e.g. MetatomicPotential).

Reimplemented from eonc::Potential.

Definition at line 73 of file LAMMPSPot.h.

                                                                     {
    return true;
  }

noteScreenGeometry()

bool LAMMPSPot::noteScreenGeometry ( long N, const double * box )

private

Definition at line 395 of file LAMMPSPot.cpp.

                                                            {
  bool changed = !screenHaveGeom_ || N != screenAtoms_;
  if (!changed) {
    for (int i = 0; i < 9; ++i) {
      if (screenBox_[i] != box[i]) {
        changed = true;
        break;
      }
    }
  }
  screenHaveGeom_ = true;
  screenAtoms_ = N;
  std::memcpy(screenBox_, box, sizeof(screenBox_));
  return changed;
}

runWorkerLoop()

void LAMMPSPot::runWorkerLoop ( )

private

Definition at line 295 of file LAMMPSPot.cpp.

                              {
  workerChild_ = true;
  // Running in the forked child.  Evaluate forces with an in-process LAMMPS
  // (this child's own MPI_COMM_WORLD) and stream results back to the parent.
  for (;;) {
    long N = 0;
    if (!readExact(reqFd, &N, sizeof(N))) {
      _exit(0); // request pipe closed -> shut down cleanly
    }
    if (N < 0) {
      _exit(0); // explicit shutdown sentinel from stopWorker()
    }
    std::vector<int> atomicNrs(static_cast<size_t>(N));
    std::vector<double> R(static_cast<size_t>(3 * N));
    double box[9];
    if (!readExact(reqFd, atomicNrs.data(),
                   sizeof(int) * static_cast<size_t>(N)) ||
        !readExact(reqFd, box, sizeof(box)) ||
        !readExact(reqFd, R.data(),
                   sizeof(double) * static_cast<size_t>(3 * N))) {
      _exit(1);
    }
    std::vector<double> mask(static_cast<size_t>(3 * N), 0.0);
    if (!readExact(reqFd, mask.data(),
                   sizeof(double) * static_cast<size_t>(3 * N))) {
      _exit(1);
    }
    setFixedMask(N, mask.data());
 
    std::vector<double> F(static_cast<size_t>(3 * N), 0.0);
    double U = 0.0;
    int status = 0;
    try {
      forceLocal(N, R.data(), atomicNrs.data(), F.data(), &U, box);
    } catch (...) {
      status = 1;
      haveStress_ = false;
    }
    double stressRaw[9] = {};
    if (haveStress_) {
      for (int i = 0; i < 9; ++i) {
        stressRaw[i] = stress_(i / 3, i % 3);
      }
    }
 
    if (!writeExact(resFd, &status, sizeof(status)) ||
        !writeExact(resFd, &U, sizeof(U)) ||
        !writeExact(resFd, F.data(),
                    sizeof(double) * static_cast<size_t>(3 * N)) ||
        !writeExact(resFd, stressRaw, sizeof(stressRaw))) {
      _exit(1);
    }
  }
}

setFixedMask()

void LAMMPSPot::setFixedMask ( long nAtoms, const double * isFixed )

overridevirtual

Optional frozen-atom mask (nAtoms*3, 1.0 = fixed).

Default no-op. Matter calls this immediately before force() so wrappers that need it (LAMMPS setforce) see Atom.fixed without changing force().

Reimplemented from eonc::Potential.

Definition at line 87 of file LAMMPSPot.cpp.

                                                               {
  std::unique_lock<std::mutex> lock(maskMutex_, std::defer_lock);
  if (!workerChild_) {
    lock.lock();
  }
  if (nAtoms <= 0 || isFixed == nullptr) {
    fixedMask_.clear();
    maskN_ = 0;
    return;
  }
  fixedMask_.assign(isFixed, isFixed + 3 * nAtoms);
  maskN_ = nAtoms;
}

stopWorker()

void LAMMPSPot::stopWorker ( )

private

Definition at line 350 of file LAMMPSPot.cpp.

                           {
  if (!workerSpawned)
    return;
  if (reqFd >= 0) {
    // Send an explicit shutdown sentinel, then close.  A sentinel (rather than
    // relying on pipe EOF) guarantees the child exits even when sibling worker
    // processes hold an inherited copy of this write end.
    long sentinel = -1;
    writeExact(reqFd, &sentinel, sizeof(sentinel));
    close(reqFd);
    reqFd = -1;
  }
  if (resFd >= 0) {
    close(resFd);
    resFd = -1;
  }
  if (workerPid > 0) {
    // A wedged worker never acts on the sentinel or the closed pipe, and an
    // unconditional wait then blocks the client forever at no CPU cost -- the
    // hang this teardown exists to avoid. Give the child a brief chance to
    // exit on its own, then insist.
    int st = 0;
    bool reaped = false;
    for (int i = 0; i < 100; ++i) { // up to ~1 s
      pid_t r = waitpid(workerPid, &st, WNOHANG);
      if (eonc::lammpsWorkerReaped(r, workerPid, errno)) {
        reaped = true;
        break;
      }
      usleep(10000);
    }
    if (!reaped) {
      kill(workerPid, SIGKILL);
      for (;;) {
        pid_t r = waitpid(workerPid, &st, 0);
        if (r == workerPid || (r < 0 && errno != EINTR)) {
          break;
        }
      }
    }
    workerPid = -1;
  }
  workerSpawned = false;
}

Member Data Documentation

fixedMask_

std::vector<double> LAMMPSPot::fixedMask_

private

Definition at line 121 of file LAMMPSPot.h.

haveStress_

bool LAMMPSPot::haveStress_ {false}

private

Definition at line 120 of file LAMMPSPot.h.

{false};

lammpsLogging_

bool LAMMPSPot::lammpsLogging_ {false}

private

Definition at line 97 of file LAMMPSPot.h.

{false};

lammpsLogIndex_

int LAMMPSPot::lammpsLogIndex_ {0}

private

Definition at line 98 of file LAMMPSPot.h.

{0};

LAMMPSObj

void* LAMMPSPot::LAMMPSObj {nullptr}

private

Definition at line 114 of file LAMMPSPot.h.

{nullptr};

lammpsScreenPath_

std::string LAMMPSPot::lammpsScreenPath_

private

Definition at line 103 of file LAMMPSPot.h.

lammpsScreenPos_

std::int64_t LAMMPSPot::lammpsScreenPos_ {0}

private

Definition at line 104 of file LAMMPSPot.h.

{0};

lammpsScreenRestart_

bool LAMMPSPot::lammpsScreenRestart_ {false}

private

Definition at line 107 of file LAMMPSPot.h.

{false};

lammpsThr

int LAMMPSPot::lammpsThr {0}

private

Definition at line 96 of file LAMMPSPot.h.

{0};

loader_

eonc::ILammpsLoader& LAMMPSPot::loader_

private

Definition at line 95 of file LAMMPSPot.h.

maskMutex_

std::mutex LAMMPSPot::maskMutex_

private

Definition at line 99 of file LAMMPSPot.h.

maskN_

long LAMMPSPot::maskN_ {0}

private

Definition at line 122 of file LAMMPSPot.h.

{0};

numberOfAtoms

long LAMMPSPot::numberOfAtoms {0}

private

Definition at line 112 of file LAMMPSPot.h.

{0};

oldBox

double LAMMPSPot::oldBox[9] {}

private

Definition at line 113 of file LAMMPSPot.h.

{};

realunits

bool LAMMPSPot::realunits {false}

private

Definition at line 118 of file LAMMPSPot.h.

{false};

reqFd

int LAMMPSPot::reqFd {-1}

private

Definition at line 154 of file LAMMPSPot.h.

{-1}; // parent writes requests here (child stdin side)

resFd

int LAMMPSPot::resFd {-1}

private

Definition at line 155 of file LAMMPSPot.h.

{-1}; // parent reads results here (child stdout side)

screenAtoms_

long LAMMPSPot::screenAtoms_ {0}

private

Definition at line 160 of file LAMMPSPot.h.

{0};

screenBox_

double LAMMPSPot::screenBox_[9] {}

private

Definition at line 161 of file LAMMPSPot.h.

{};

screenHaveGeom_

bool LAMMPSPot::screenHaveGeom_ {false}

private

Definition at line 159 of file LAMMPSPot.h.

{false};

stress_

Matrix3d LAMMPSPot::stress_ {Matrix3d::Zero()}

private

Definition at line 119 of file LAMMPSPot.h.

{Matrix3d::Zero()};

workerChild_

bool LAMMPSPot::workerChild_ {false}

private

Definition at line 108 of file LAMMPSPot.h.

{false};

workerMutex

std::mutex LAMMPSPot::workerMutex

private

Definition at line 126 of file LAMMPSPot.h.

workerPid

int LAMMPSPot::workerPid {-1}

private

Definition at line 153 of file LAMMPSPot.h.

{-1};

workerRespawnsLeft

int LAMMPSPot::workerRespawnsLeft {3}

private

Definition at line 144 of file LAMMPSPot.h.

{3};

workerSpawned

bool LAMMPSPot::workerSpawned {false}

private

Definition at line 156 of file LAMMPSPot.h.

{false};

---

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

• /home/runner/work/eOn/eOn/include/eon/potentials/LAMMPS/LAMMPSPot.h
• /home/runner/work/eOn/eOn/client/potentials/LAMMPS/LAMMPSPot.cpp