#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 Parameters &p)
	~LAMMPSPot ()
void	cleanMemory ()
void	force (long N, const double R, const int atomicNrs, double F, double U, double variance, const double box) 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	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	makeNewLAMMPS (long N, const double R, const int atomicNrs, const double *box)
void	ensureWorker ()
void	runWorkerLoop ()
void	stopWorker ()
void	forceLocal (long N, const double R, const int atomicNrs, double F, double U, const double *box)

Private Attributes
int	lammpsThr {0}
long	numberOfAtoms {0}
double	oldBox [9] {}
void *	LAMMPSObj {nullptr}
bool	realunits {false}
int	workerRespawnsLeft {3}
std::mutex	workerMutex
int	workerPid {-1}
int	reqFd {-1}
int	resFd {-1}
bool	workerSpawned {false}

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

Detailed Description

Definition at line 21 of file LAMMPSPot.h.

Constructor & Destructor Documentation

◆ LAMMPSPot()

LAMMPSPot::LAMMPSPot ( const Parameters & p )

Definition at line 40 of file LAMMPSPot.cpp.

    : Potential(p),
      lammpsThr{p.potential_options.LAMMPSThreads}
#ifdef EONMPI
      ,
      mpiComm{p.potential_options.MPIClientComm}
#endif
{
  // Fail fast if LAMMPS library not available
  eonc::LammpsLoader::instance().require_loaded();
#if !defined(EONMPI) && !defined(IS_WINDOWS)
  // 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();
#endif
}

◆ ~LAMMPSPot()

LAMMPSPot::~LAMMPSPot ( )

Definition at line 61 of file LAMMPSPot.cpp.

61{ cleanMemory(); }

LAMMPSPot::cleanMemory

void cleanMemory()

Definition LAMMPSPot.cpp:63

Member Function Documentation

◆ cleanMemory()

void LAMMPSPot::cleanMemory ( )

Definition at line 63 of file LAMMPSPot.cpp.

                            {
#if !defined(EONMPI) && !defined(IS_WINDOWS)
  stopWorker();
#endif
  if (LAMMPSObj != nullptr) {
    eonc::LammpsLoader::instance().close(LAMMPSObj);
    LAMMPSObj = nullptr;
  }
}

◆ ensureWorker()

void LAMMPSPot::ensureWorker ( )

private

Definition at line 140 of file LAMMPSPot.cpp.

                             {
  if (workerSpawned)
    return;
 
  int reqPipe[2]; // parent -> child
  int resPipe[2]; // child -> parent
  if (pipe(reqPipe) != 0 || pipe(resPipe) != 0) {
    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) {
    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 273 of file LAMMPSPot.cpp.

                                                                      {
  variance = nullptr;
 
#ifdef EONMPI
  forceLocal(N, R, atomicNrs, F, U, box);
#elif defined(IS_WINDOWS)
  // No fork/pipe on Windows; call forceLocal directly.
  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;
  }
  ensureWorker();
 
  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))) {
    // 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;
  if (!readExact(resFd, &status, sizeof(status)) ||
      !readExact(resFd, U, sizeof(double)) ||
      !readExact(resFd, F, sizeof(double) * static_cast<size_t>(3 * N))) {
    --workerRespawnsLeft;
    EONC_LOG_WARNING(
        "[LAMMPSPot] worker died during force eval; {} respawns left "
        "(eon-7416)",
        workerRespawnsLeft);
    stopWorker();
    rejectGeometry(U, F, N);
    return;
  }
  if (status != 0) {
    --workerRespawnsLeft;
    EONC_LOG_WARNING(
        "[LAMMPSPot] worker reported an evaluation error; {} respawns left "
        "(eon-7416)",
        workerRespawnsLeft);
    stopWorker();
    rejectGeometry(U, F, N);
    return;
  }
  // 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 388 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 = eonc::LammpsLoader::instance();
 
    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));
    lmp.command(LAMMPSObj, "run 1 pre no post no");
 
    auto *pe =
        static_cast<double *>(lmp.extract_variable(LAMMPSObj, "pe", nullptr));
    *U = *pe;
    free(pe);
 
    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"));
 
    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);
  } catch (...) {
    fpeh.restore_fpe();
    throw;
  }
  fpeh.restore_fpe();
}

◆ makeNewLAMMPS()

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

private

Definition at line 453 of file LAMMPSPot.cpp.

                                                 {
  auto &lmp = eonc::LammpsLoader::instance();
 
  numberOfAtoms = N;
  std::memcpy(oldBox, box, 9 * sizeof(double));
 
  if (LAMMPSObj != nullptr) {
    eonc::LammpsLoader::instance().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 char *lmpargv[] = {"liblammps", "-log", "none",    "-echo", "log",
                           "-screen",   "none", "-suffix", "omp"};
  int lmpargc = sizeof(lmpargv) / sizeof(const char *);
  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, const_cast<char **>(lmpargv), inst_comm, nullptr);
#else
  const char *lmpargv[] = {"liblammps", "-log",    "none", "-echo",
                           "log",       "-screen", "none"};
  int lmpargc = sizeof(lmpargv) / sizeof(const char *);
  LAMMPSObj = lmp.open_no_mpi(lmpargc, const_cast<char **>(lmpargv), nullptr);
#endif
 
  if (lammpsThr > 0) {
    std::string cmd = std::format("package omp {} force/neigh", lammpsThr);
    lmp.command(LAMMPSObj, 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 {
    EONC_LOG_ERROR("[LAMMPS] in.lammps not found in working directory");
    return;
  }
 
  if (realunits) {
    lmp.command(LAMMPSObj, "units real");
  } else {
    lmp.command(LAMMPSObj, "units metal");
  }
 
  lmp.command(LAMMPSObj, "atom_style charge");
  lmp.command(LAMMPSObj, "atom_modify map array sort 0 0");
  lmp.command(LAMMPSObj, "neigh_modify delay 1");
 
  // Define periodic cell (prism for non-orthorhombic)
  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]);
  lmp.command(LAMMPSObj, region_cmd.c_str());
 
  std::string create_box_cmd = std::format("create_box {} cell", ntypes);
  lmp.command(LAMMPSObj, 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);
    lmp.command(LAMMPSObj, atom_cmd.c_str());
  }
 
  lmp.command(LAMMPSObj, "mass * 1.0");
 
  // Load user LAMMPS input script
  lmp.file(LAMMPSObj, "in.lammps");
 
  // Define variables for force/energy extraction
  lmp.command(LAMMPSObj, "variable fx atom fx");
  lmp.command(LAMMPSObj, "variable fy atom fy");
  lmp.command(LAMMPSObj, "variable fz atom fz");
  lmp.command(LAMMPSObj, "variable pe equal pe");
}

◆ 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 24 of file LAMMPSPot.h.

                                                                     {
    return true;
  }

◆ runWorkerLoop()

void LAMMPSPot::runWorkerLoop ( )

private

Definition at line 192 of file LAMMPSPot.cpp.

                              {
  // 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> 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;
    }
 
    if (!writeExact(resFd, &status, sizeof(status)) ||
        !writeExact(resFd, &U, sizeof(U)) ||
        !writeExact(resFd, F.data(),
                    sizeof(double) * static_cast<size_t>(3 * N))) {
      _exit(1);
    }
  }
}

◆ stopWorker()

void LAMMPSPot::stopWorker ( )

private

Definition at line 232 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 (r == workerPid || r < 0) {
        reaped = true;
        break;
      }
      usleep(10000);
    }
    if (!reaped) {
      kill(workerPid, SIGKILL);
      waitpid(workerPid, &st, 0);
    }
    workerPid = -1;
  }
  workerSpawned = false;
}

Member Data Documentation

◆ LAMMPSObj

void* LAMMPSPot::LAMMPSObj {nullptr}

private

Definition at line 40 of file LAMMPSPot.h.

40{nullptr};

◆ lammpsThr

int LAMMPSPot::lammpsThr {0}

private

Definition at line 34 of file LAMMPSPot.h.

34{0};

◆ numberOfAtoms

long LAMMPSPot::numberOfAtoms {0}

private

Definition at line 38 of file LAMMPSPot.h.

38{0};

◆ oldBox

double LAMMPSPot::oldBox[9] {}

private

Definition at line 39 of file LAMMPSPot.h.

39{};

◆ realunits

bool LAMMPSPot::realunits {false}

private

Definition at line 43 of file LAMMPSPot.h.

43{false};

◆ reqFd

int LAMMPSPot::reqFd {-1}

private

Definition at line 72 of file LAMMPSPot.h.

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

◆ resFd

int LAMMPSPot::resFd {-1}

private

Definition at line 73 of file LAMMPSPot.h.

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

◆ workerMutex

std::mutex LAMMPSPot::workerMutex

private

Definition at line 70 of file LAMMPSPot.h.

◆ workerPid

int LAMMPSPot::workerPid {-1}

private

Definition at line 71 of file LAMMPSPot.h.

71{-1};

◆ workerRespawnsLeft

int LAMMPSPot::workerRespawnsLeft {3}

private

Definition at line 61 of file LAMMPSPot.h.

61{3};

◆ workerSpawned

bool LAMMPSPot::workerSpawned {false}

private

Definition at line 74 of file LAMMPSPot.h.

74{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

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ LAMMPSPot()

◆ ~LAMMPSPot()

Member Function Documentation

◆ cleanMemory()

◆ ensureWorker()

◆ force()

◆ forceLocal()

◆ makeNewLAMMPS()

◆ needsPerImageInstance()

◆ runWorkerLoop()

◆ stopWorker()

Member Data Documentation

◆ LAMMPSObj

◆ lammpsThr

◆ numberOfAtoms

◆ oldBox

◆ realunits

◆ reqFd

◆ resFd

◆ workerMutex

◆ workerPid

◆ workerRespawnsLeft

◆ workerSpawned