ConFileIO.cpp

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Copyright (c) 2010--present, eOn Development Team
** All rights reserved.
**
** Repo:
** https://github.com/TheochemUI/eOn
*/
#include "eon/ConFileIO.h"
#include "ReadconDbMirror.h"
#include "eon/Eigen.h"
#include "eon/EonLogger.h"
#include "eon/HelperFunctions.h"
#include "eon/Matter.h"
#include "eon/SafeMath.h"
 
#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <filesystem>
#include <format>
#include <fstream>
#include <mutex>
#include <numeric>
#include <optional>
#include <stdexcept>
#include <string>
#include <system_error>
#include <unordered_map>
#include <vector>
 
namespace {
 
using Matter = eonc::Matter;
namespace fs = std::filesystem;
 
constexpr uint8_t kConPrecision = 17;
constexpr uint8_t kConvelPrecision = 6;
 
std::string strip_nl(const std::string &s) {
  std::string str(s);
  while (!str.empty() && (str.back() == '\n' || str.back() == '\r'))
    str.pop_back();
  return str;
}
 
std::string canonical_generator_header(const std::string &header) {
  auto stripped = strip_nl(header);
  if (!stripped.empty()) {
    return stripped;
  }
  return "Generated by eOn";
}
 
std::string ensure_extension(std::string filename, std::string_view ext) {
  fs::path path(filename);
  const auto name = path.filename().string();
  const bool has_compound =
      name.size() > ext.size() && (name.ends_with(std::string(ext) + ".gz") ||
                                   name.ends_with(std::string(ext) + ".zst"));
  if (!has_compound && path.extension() != ext) {
    path += ext;
  }
  return path.string();
}
 
std::string symbol_for_z(long atomic_nr) {
  return readcon::z_to_symbol(static_cast<uint64_t>(atomic_nr));
}
 
std::vector<double> flat_row_major(const AtomMatrix &m) {
  static_assert(AtomMatrix::IsRowMajor,
                "flat_row_major assumes AtomMatrix RowMajor Nx3 layout");
  const auto n = static_cast<size_t>(m.rows());
  std::vector<double> flat(n * 3);
  std::memcpy(flat.data(), m.data(), flat.size() * sizeof(double));
  return flat;
}
 
void apply_frame_metadata(readcon::ConFrameBuilder &builder,
                          const eonc::io::ConFrameMetadata *metadata) {
  if (metadata == nullptr) {
    return;
  }
  if (metadata->raw_json && !metadata->raw_json->empty()) {
    builder.set_metadata_json(*metadata->raw_json);
  }
  if (metadata->energy) {
    builder.set_energy(*metadata->energy);
  }
  if (metadata->frame_index) {
    builder.set_frame_index(*metadata->frame_index);
  }
  if (metadata->time) {
    builder.set_time(*metadata->time);
  }
  if (metadata->timestep) {
    builder.set_timestep(*metadata->timestep);
  }
  if (metadata->neb_bead) {
    builder.set_neb_bead(*metadata->neb_bead);
  }
  if (metadata->neb_band) {
    builder.set_neb_band(*metadata->neb_band);
  }
  if (metadata->potential_type && !metadata->potential_type->empty()) {
    builder.set_string_metadata("potential_type", *metadata->potential_type);
  }
  for (const auto &[key, value] : metadata->scalars) {
    builder.set_scalar_metadata(key, value);
  }
  for (const auto &[key, value] : metadata->strings) {
    builder.set_string_metadata(key, value);
  }
}
 
eonc::io::IoStatus write_frames(const fs::path &path,
                                const std::vector<readcon::ConFrame> &frames,
                                uint8_t precision) {
  try {
    const auto compression =
        readcon::ConFrameWriter::compression_from_extension(path);
    readcon::ConFrameWriter writer(path, compression, precision);
    writer.extend(frames);
    return eonc::io::IoStatus::Ok;
  } catch (const std::exception &e) {
    EONC_LOG_ERROR("Failed to write {}: {}", path.string(), e.what());
    return eonc::io::IoStatus::WriteError;
  }
}

std::vector<size_t> matter_order(const std::vector<readcon::Atom> &atoms) {
  const size_t n = atoms.size();
  std::vector<size_t> order(n);
  std::iota(order.begin(), order.end(), size_t{0});
  if (n < 2) {
    return order;
  }
  const bool ascending =
      std::is_sorted(atoms.begin(), atoms.end(),
                     [](const readcon::Atom &a, const readcon::Atom &b) {
                       return a.atom_id < b.atom_id;
                     }) &&
      std::adjacent_find(atoms.begin(), atoms.end(),
                         [](const readcon::Atom &a, const readcon::Atom &b) {
                           return a.atom_id == b.atom_id;
                         }) == atoms.end();
  if (ascending) {
    return order;
  }
  std::vector<uint64_t> ids;
  ids.reserve(n);
  for (const auto &atom : atoms) {
    ids.push_back(atom.atom_id);
  }
  std::vector<uint64_t> sorted(ids);
  std::sort(sorted.begin(), sorted.end());
  if (std::adjacent_find(sorted.begin(), sorted.end()) != sorted.end()) {
    return order;
  }
  std::stable_sort(order.begin(), order.end(),
                   [&ids](size_t a, size_t b) { return ids[a] < ids[b]; });
  return order;
}

struct FileStamp {
  std::uintmax_t size{0};
  fs::file_time_type mtime{};
  friend bool operator==(const FileStamp &, const FileStamp &) = default;
};
 
std::optional<FileStamp> stamp_file(const fs::path &path) {
  std::error_code ec;
  const auto size = fs::file_size(path, ec);
  if (ec) {
    return std::nullopt;
  }
  const auto mtime = fs::last_write_time(path, ec);
  if (ec) {
    return std::nullopt;
  }
  return FileStamp{size, mtime};
}

std::string append_key(const fs::path &path) {
  std::error_code ec;
  auto resolved = fs::weakly_canonical(path, ec);
  if (ec || resolved.empty()) {
    resolved = fs::absolute(path, ec);
    if (ec) {
      return path.lexically_normal().string();
    }
  }
  return resolved.lexically_normal().string();
}

std::mutex &append_mutex() {
  static std::mutex m;
  return m;
}
 
std::unordered_map<std::string, FileStamp> &append_stamps() {
  static std::unordered_map<std::string, FileStamp> stamps;
  return stamps;
}

void remember_stamp(const std::string &key, const fs::path &path) {
  if (auto stamp = stamp_file(path)) {
    append_stamps()[key] = *stamp;
  } else {
    append_stamps().erase(key);
  }
}

bool tail_is_ours(const std::string &key, const fs::path &path) {
  const auto it = append_stamps().find(key);
  if (it == append_stamps().end()) {
    return false;
  }
  const auto stamp = stamp_file(path);
  return stamp.has_value() && *stamp == it->second;
}

std::optional<long> last_xyz_natoms(const fs::path &path) {
  std::ifstream in(path);
  if (!in) {
    return std::nullopt;
  }
  std::optional<long> last;
  std::string line;
  while (std::getline(in, line)) {
    if (line.empty()) {
      continue;
    }
    long n = 0;
    try {
      n = std::stol(line);
    } catch (const std::exception &) {
      return std::nullopt;
    }
    if (n < 0) {
      return std::nullopt;
    }
    if (!std::getline(in, line)) {
      return std::nullopt;
    }
    for (long i = 0; i < n; ++i) {
      if (!std::getline(in, line)) {
        return std::nullopt;
      }
    }
    last = n;
  }
  return last;
}

bool ends_with_newline(const fs::path &path) {
  std::ifstream in(path, std::ios::binary | std::ios::ate);
  if (!in) {
    return true;
  }
  const auto len = static_cast<std::streamoff>(in.tellg());
  if (len <= 0) {
    return true;
  }
  in.seekg(len - 1);
  char last = '\n';
  if (!in.get(last)) {
    return true;
  }
  return last == '\n';
}

eonc::io::IoStatus append_frames(const fs::path &path,
                                 const std::vector<readcon::ConFrame> &frames,
                                 uint8_t precision) {
  static std::atomic<uint64_t> scratch_counter{0};
  const auto scratch =
      path.parent_path() /
      std::format(".{}.eon-append-{}.tmp", path.filename().string(),
                  scratch_counter.fetch_add(1, std::memory_order_relaxed));
 
  std::string bytes;
  try {
    {
      readcon::ConFrameWriter writer(
          scratch, readcon::ConFrameWriter::Compression::None, precision);
      writer.extend(frames);
    }
    std::ifstream in(scratch, std::ios::binary | std::ios::ate);
    if (!in) {
      throw std::runtime_error("serialized frame not readable back");
    }
    const auto len = static_cast<std::streamoff>(in.tellg());
    if (len <= 0) {
      throw std::runtime_error("serialized frame is empty");
    }
    bytes.resize(static_cast<size_t>(len));
    in.seekg(0);
    if (!in.read(bytes.data(), static_cast<std::streamsize>(len))) {
      throw std::runtime_error("short read of serialized frame");
    }
  } catch (const std::exception &e) {
    std::error_code ec;
    fs::remove(scratch, ec);
    EONC_LOG_ERROR("Failed to serialize frame for {}: {}", path.string(),
                   e.what());
    return eonc::io::IoStatus::WriteError;
  }
  std::error_code ec;
  fs::remove(scratch, ec);
 
  std::ofstream out(path, std::ios::binary | std::ios::app);
  if (!out) {
    EONC_LOG_ERROR("Failed to open {} for append", path.string());
    return eonc::io::IoStatus::WriteError;
  }
  // A frame header must start its own line. eOn's own frames end in a
  // newline; a hand-written or foreign tail may not.
  if (!ends_with_newline(path)) {
    out.put('\n');
  }
  out.write(bytes.data(), static_cast<std::streamsize>(bytes.size()));
  out.flush();
  if (!out) {
    EONC_LOG_ERROR("Failed to append to {}", path.string());
    return eonc::io::IoStatus::WriteError;
  }
  return eonc::io::IoStatus::Ok;
}

readcon::ConFrameBuilder seed_builder(eonc::Matter &m,
                                      const std::array<std::string, 2> &prebox,
                                      const std::array<std::string, 2> &postbox,
                                      const std::vector<uint64_t> &atom_ids) {
  auto [lengths, angles_deg] = eonc::io::cell_to_lengths_angles(m);
  readcon::ConFrameBuilder builder(
      {lengths[0], lengths[1], lengths[2]},
      {angles_deg[0], angles_deg[1], angles_deg[2]}, prebox, postbox);
 
  const long n = m.numberOfAtoms();
  if (atom_ids.size() != static_cast<size_t>(n)) {
    throw std::invalid_argument("seed_builder: atom_ids size mismatch");
  }
  for (long i = 0; i < n; ++i) {
    const auto mask = m.getFixedMask(i);
    builder.add_atom(symbol_for_z(m.getAtomicNr(i)), 0.0, 0.0, 0.0, mask,
                     atom_ids[static_cast<size_t>(i)], m.getMass(i));
  }
  return builder;
}
 
bool should_write_forces(const Matter &m,
                         const eonc::io::ConFrameMetadata *metadata) {
  // Per-write metadata wins so two callers can disagree without a race on
  // the process-wide flag. Parameters next, so pyeonclient setting the
  // field on the bound Parameters object is enough. The atomic global is
  // last (INI parse and the explicit setter).
  if (metadata && metadata->write_con_forces.has_value()) {
    return *metadata->write_con_forces;
  }
  if (m.getWriteConForces()) {
    return true;
  }
  return eonc::io::write_con_forces();
}
 
void apply_geometry(readcon::ConFrameBuilder &builder, Matter &m,
                    bool with_velocities,
                    const eonc::io::ConFrameMetadata *metadata) {
  const long n = m.numberOfAtoms();
  if (n <= 0) {
    return;
  }
 
  // Prefer bulk flat setters (declared sections, no pointer-lifetime hazards
  // under ASAN). AtomMatrix is RowMajor Nx3 matching the flat layout.
  builder.set_positions_from_flat(flat_row_major(m.getPositions()));
 
  // getForcesRaw() may trigger pot evaluation if recomputePotential is set;
  // only enter when the pot is already clean so we serialize cached forces.
  // Gated: force sections are opt-in ([Main] write_con_forces) because
  // ASE-class readers reject frames that carry them.
  if (should_write_forces(m, metadata) && !m.needsForceUpdate()) {
    builder.set_forces_from_flat(flat_row_major(m.getForcesRaw()));
  }
 
  if (metadata != nullptr && !metadata->displacements.empty()) {
    if (metadata->displacements.size() != static_cast<size_t>(3 * n)) {
      throw std::invalid_argument("displacements size is not 3 x atoms");
    }
    builder.set_displacements_from_flat(metadata->displacements);
  }
  if (metadata != nullptr && !metadata->spreads.empty()) {
    if (metadata->spreads.size() != static_cast<size_t>(3 * n)) {
      throw std::invalid_argument("spreads size is not 3 x atoms");
    }
    builder.set_spreads_from_flat(metadata->spreads);
  }
 
  if (with_velocities) {
    const AtomMatrix vel = m.getVelocities();
    for (long i = 0; i < n; ++i) {
      builder.set_atom_velocity(static_cast<size_t>(i),
                                {vel(i, 0), vel(i, 1), vel(i, 2)});
    }
  }
}
 
void collect_ids_headers(Matter &m, std::vector<uint64_t> &atom_ids,
                         std::array<std::string, 2> &prebox,
                         std::array<std::string, 2> &postbox) {
  const long n = m.numberOfAtoms();
  atom_ids.resize(static_cast<size_t>(n));
  for (long i = 0; i < n; ++i) {
    atom_ids[static_cast<size_t>(i)] = static_cast<uint64_t>(m.getAtomIndex(i));
  }
  const auto &hdr = m.getHeaderCon();
  prebox = {canonical_generator_header(hdr[0]), strip_nl(hdr[1])};
  postbox = {strip_nl(hdr[3]), strip_nl(hdr[4])};
}
 
readcon::ConFrame frame_from_matter(Matter &m,
                                    const eonc::io::ConFrameMetadata *metadata,
                                    bool with_velocities) {
  std::vector<uint64_t> atom_ids;
  std::array<std::string, 2> prebox;
  std::array<std::string, 2> postbox;
  collect_ids_headers(m, atom_ids, prebox, postbox);
 
  auto builder = seed_builder(m, prebox, postbox, atom_ids);
 
  eonc::io::ConFrameMetadata auto_meta;
  const eonc::io::ConFrameMetadata *meta_ptr = metadata;
  if (!m.needsForceUpdate()) {
    if (metadata == nullptr) {
      auto_meta.energy = m.getPotentialEnergy();
      meta_ptr = &auto_meta;
    } else if (!metadata->energy) {
      auto_meta = *metadata;
      auto_meta.energy = m.getPotentialEnergy();
      meta_ptr = &auto_meta;
    }
  }
  apply_frame_metadata(builder, meta_ptr);
  apply_geometry(builder, m, with_velocities, meta_ptr);
  return builder.build();
}
 
} // namespace
 
namespace eonc::io {
 
namespace {
// Process-wide opt-in for force sections in written frames; classic con
// layout (no sections) keeps ASE-class readers working on our outputs.
std::atomic<bool> g_write_con_forces{false};
} // namespace
 
readcon::ConFrame matterToConFrame(Matter &m,
                                   const ConFrameMetadata *metadata) {
  return frame_from_matter(m, metadata, /*with_velocities=*/false);
}
 
void set_write_con_forces(bool enabled) noexcept {
  g_write_con_forces.store(enabled, std::memory_order_relaxed);
}
 
bool write_con_forces() noexcept {
  return g_write_con_forces.load(std::memory_order_relaxed);
}
 
ConFrameMetadata metadata_from_frame(const readcon::ConFrame &frame) {
  ConFrameMetadata meta;
  meta.energy = frame.energy_opt();
  meta.frame_index = frame.frame_index_opt();
  meta.time = frame.time_opt();
  meta.timestep = frame.timestep_opt();
  meta.neb_bead = frame.neb_bead_opt();
  meta.neb_band = frame.neb_band_opt();
  meta.potential_type = frame.potential_type();
  const auto json = frame.metadata_json();
  if (!json.empty() && json != "{}") {
    meta.raw_json = json;
  }
  return meta;
}
 
std::pair<std::array<double, 3>, std::array<double, 3>>
cell_to_lengths_angles(const Matter &m) {
  const Matrix3d cell = m.getCell();
  std::array<double, 3> lengths;
  lengths[0] = cell.row(0).norm();
  lengths[1] = cell.row(1).norm();
  lengths[2] = cell.row(2).norm();
  // CON header line 4 is alpha beta gamma:
  // alpha = angle(b,c), beta = angle(a,c), gamma = angle(a,b).
  std::array<double, 3> angles;
  angles[0] = eonc::safemath::safe_acos(eonc::safemath::safe_div(
                  cell.row(1).dot(cell.row(2)), lengths[1] * lengths[2])) *
              180.0 / eonc::helpers::pi;
  angles[1] = eonc::safemath::safe_acos(eonc::safemath::safe_div(
                  cell.row(0).dot(cell.row(2)), lengths[0] * lengths[2])) *
              180.0 / eonc::helpers::pi;
  angles[2] = eonc::safemath::safe_acos(eonc::safemath::safe_div(
                  cell.row(0).dot(cell.row(1)), lengths[0] * lengths[1])) *
              180.0 / eonc::helpers::pi;
  return {lengths, angles};
}
 
void resetConAppendState() {
  const std::lock_guard<std::mutex> guard(append_mutex());
  append_stamps().clear();
}
 
IoStatus matter2con(Matter &m, std::string filename, bool append,
                    const ConFrameMetadata *metadata) {
  filename = ensure_extension(std::move(filename), ".con");
 
  m.applyPeriodicBoundaryIfEnabled();
 
  const fs::path path(filename);
  const auto compression =
      readcon::ConFrameWriter::compression_from_extension(path);
  const bool streamable =
      compression == readcon::ConFrameWriter::Compression::None;
 
  // The lock covers the stamp table and the read-then-extend window on one
  // path. Nothing under it evaluates a potential: frame_from_matter reads
  // cached energy and forces only while the pot is clean.
  const std::lock_guard<std::mutex> guard(append_mutex());
  const auto key = append_key(path);
  const bool exists = fs::exists(path);
  const bool concatenate = append && exists && streamable;
 
  // ConFrame is move-only (no default ctor), so frames carries the new frame
  // either way. A gzip or zstd target cannot be concatenated: readcon-core
  // reads a single gzip member, so appended members would be invisible
  // on read-back, and the streaming writers cannot be flushed frame by frame.
  // Those targets keep the whole-file rewrite, which stays O(N) per append.
  std::vector<readcon::ConFrame> frames;
  if (append && exists && !streamable) {
    try {
      // Prefer read_all_frames (single ownership hand-off) over the iterator
      // path for append rewrite; simpler lifetime for ASAN/CI envs.
      frames = readcon::read_all_frames(path);
    } catch (const std::exception &e) {
      EONC_LOG_ERROR("Failed to append to {}: {}", filename, e.what());
      return IoStatus::AppendError;
    }
  } else if (concatenate && !tail_is_ours(key, path)) {
    // Refuse to extend a file eOn cannot parse, matching the rewrite path:
    // the target keeps its bytes and the caller sees AppendError. Frames this
    // process wrote and nobody touched since need no such check.
    try {
      [[maybe_unused]] const auto existing = readcon::read_all_frames(path);
    } catch (const std::exception &e) {
      EONC_LOG_ERROR("Failed to append to {}: {}", filename, e.what());
      return IoStatus::AppendError;
    }
  }
  try {
    frames.push_back(frame_from_matter(m, metadata, /*with_velocities=*/false));
  } catch (const std::exception &e) {
    EONC_LOG_ERROR("Failed to build frame for {}: {}", filename, e.what());
    return IoStatus::InvalidArgument;
  }
 
  const auto status = concatenate ? append_frames(path, frames, kConPrecision)
                                  : write_frames(path, frames, kConPrecision);
  if (io_ok(status)) {
    remember_stamp(key, path);
    if (!concatenate) {
      mirror_con_corpus(path.string());
    }
  } else {
    append_stamps().erase(key);
  }
  return status;
}
 
IoStatus con2matter(Matter &m, std::string filename) {
  filename = ensure_extension(std::move(filename), ".con");
  try {
    auto frame = readcon::read_first_frame(filename);
    const IoStatus status = con2matter(m, frame, nullptr);
    if (io_ok(status)) {
      mirror_con_corpus(filename);
    }
    return status;
  } catch (const std::exception &e) {
    EONC_LOG_ERROR("Failed to read {}: {}", filename, e.what());
    return IoStatus::ReadError;
  }
}
 
IoStatus con2matter(Matter &m, const readcon::ConFrame &frame,
                    ConFrameMetadata *out_metadata) {
  try {
    const auto &atoms = frame.atoms();
    const auto &lengths = frame.cell();
    const auto &angles_deg = frame.angles();
    const auto &prebox = frame.prebox_header();
    const auto &postbox = frame.postbox_header();
 
    m.headerCon[0] = prebox[0] + "\n";
    m.headerCon[1] = prebox[1] + "\n";
    m.headerCon[3] = postbox[0] + "\n";
    m.headerCon[4] = postbox[1] + "\n";
 
    double angles[3] = {angles_deg[0], angles_deg[1], angles_deg[2]};
    if (angles[0] == 90.0 && angles[1] == 90.0 && angles[2] == 90.0) {
      Matrix3d cell = Matrix3d::Zero();
      cell(0, 0) = lengths[0];
      cell(1, 1) = lengths[1];
      cell(2, 2) = lengths[2];
      m.setCell(cell);
    } else {
      angles[0] *= eonc::helpers::pi / 180.0;
      angles[1] *= eonc::helpers::pi / 180.0;
      angles[2] *= eonc::helpers::pi / 180.0;
      const double alpha = angles[0];
      const double beta = angles[1];
      const double gamma = angles[2];
 
      Matrix3d cell = Matrix3d::Zero();
      cell(0, 0) = 1.0;
      cell(1, 0) = cos(gamma);
      cell(1, 1) = sin(gamma);
      cell(2, 0) = cos(beta);
      cell(2, 1) = (cos(alpha) - cell(1, 0) * cell(2, 0)) / cell(1, 1);
      cell(2, 2) = eonc::safemath::safe_sqrt(1.0 - pow(cell(2, 0), 2) -
                                             pow(cell(2, 1), 2));
 
      cell(0, 0) *= lengths[0];
      cell(1, 0) *= lengths[1];
      cell(1, 1) *= lengths[1];
      cell(2, 0) *= lengths[2];
      cell(2, 1) *= lengths[2];
      cell(2, 2) *= lengths[2];
      m.setCell(cell);
    }
    m.headerCon[2] =
        std::format("{} {} {}\n", angles_deg[0], angles_deg[1], angles_deg[2]);
 
    const auto n = static_cast<Eigen::Index>(atoms.size());
    m.resize(static_cast<long>(atoms.size()));
 
    // Undo the species grouping the .con format imposes, so an index into
    // Matter addresses the same atom as the matching row of mode.dat.
    const std::vector<size_t> order = matter_order(atoms);
    std::vector<long> file_to_matter(static_cast<size_t>(n));
    for (Eigen::Index i = 0; i < n; ++i) {
      file_to_matter[order[static_cast<size_t>(i)]] = static_cast<long>(i);
    }
    m.setFileToMatter(std::move(file_to_matter));
 
    AtomMatrix positions = AtomMatrix::Zero(n, 3);
    AtomMatrix forces = AtomMatrix::Zero(n, 3);
    AtomMatrix velocities = AtomMatrix::Zero(n, 3);
    VectorXd masses = VectorXd::Zero(n);
    VectorXi atomic_nrs = VectorXi::Zero(n);
    bool any_force = false;
    bool any_velocity = false;
 
    for (Eigen::Index i = 0; i < n; ++i) {
      const auto &atom = atoms[order[static_cast<size_t>(i)]];
      positions(i, 0) = atom.x;
      positions(i, 1) = atom.y;
      positions(i, 2) = atom.z;
      masses(i) = atom.mass;
      atomic_nrs(i) = static_cast<int>(atom.atomic_number);
      const auto fixed = atom.fixed_mask();
      m.setFixedMask(static_cast<long>(i), {fixed[0], fixed[1], fixed[2]});
      m.setAtomIndex(static_cast<long>(i),
                     static_cast<std::int64_t>(atom.atom_id));
 
      if (auto vel = atom.velocity()) {
        any_velocity = true;
        velocities(i, 0) = (*vel)[0];
        velocities(i, 1) = (*vel)[1];
        velocities(i, 2) = (*vel)[2];
      }
      if (auto force = atom.force()) {
        any_force = true;
        forces(i, 0) = (*force)[0];
        forces(i, 1) = (*force)[1];
        forces(i, 2) = (*force)[2];
      }
    }
 
    m.setMasses(masses);
    m.setAtomicNrs(atomic_nrs);
    m.setPositions(positions);
 
    if (any_velocity || frame.has_velocities()) {
      m.setVelocities(velocities);
    }
 
    const auto meta = metadata_from_frame(frame);
    if (out_metadata != nullptr) {
      *out_metadata = meta;
    }
 
    // Trust file energy+forces only when both are present. Energy-only must not
    // mark the pot clean with a zero force matrix (optimizer footgun). Prefer
    // writing raw forces (friend) so fixed-atom components survive RT; then
    // mark clean without setComputedPotential's net-force adjustment on zeros.
    const bool has_force_section = any_force || frame.has_forces();
    if (meta.energy && has_force_section) {
      m.restoreFileForces(forces, true, *meta.energy);
    } else if (has_force_section) {
      m.restoreFileForces(forces, false, 0.0);
    } else {
      // Classic geometry-only files: always recompute pot (main-era behavior).
      m.recomputePotential = true;
    }
 
    // setPositions already applied PBC when enabled; no second wrap here.
    return IoStatus::Ok;
  } catch (const std::exception &e) {
    EONC_LOG_ERROR("Failed to convert frame to matter: {}", e.what());
    return IoStatus::ReadError;
  }
}
 
IoStatus matter2convel(Matter &m, std::string filename) {
  filename = ensure_extension(std::move(filename), ".convel");
 
  m.applyPeriodicBoundaryIfEnabled();
 
  try {
    auto frame = frame_from_matter(m, nullptr, /*with_velocities=*/true);
    std::vector<readcon::ConFrame> frames;
    frames.push_back(std::move(frame));
    const IoStatus status = write_frames(filename, frames, kConvelPrecision);
    if (io_ok(status)) {
      mirror_con_corpus(filename);
    }
    return status;
  } catch (const std::exception &e) {
    EONC_LOG_ERROR("Failed to write convel {}: {}", filename, e.what());
    return IoStatus::WriteError;
  }
}
 
IoStatus convel2matter(Matter &m, std::string filename) {
  filename = ensure_extension(std::move(filename), ".convel");
  try {
    auto frame = readcon::read_first_frame(filename);
    const IoStatus status = con2matter(m, frame, nullptr);
    if (io_ok(status)) {
      mirror_con_corpus(filename);
    }
    return status;
  } catch (const std::exception &e) {
    EONC_LOG_ERROR("Failed to read convel {}: {}", filename, e.what());
    return IoStatus::ReadError;
  }
}
 
IoStatus matter2xyz(Matter &m, std::string filename, bool append) {
  filename = ensure_extension(std::move(filename), ".xyz");
 
  m.applyPeriodicBoundaryIfEnabled();
  const long n = m.numberOfAtoms();
 
  if (append && fs::exists(filename)) {
    std::error_code ec;
    const auto sz = fs::file_size(filename, ec);
    if (!ec && sz > 0) {
      const auto prev = last_xyz_natoms(filename);
      if (!prev) {
        EONC_LOG_ERROR("matter2xyz: cannot parse existing {}", filename);
        return IoStatus::AppendError;
      }
      if (*prev != n) {
        EONC_LOG_ERROR(
            "matter2xyz: append atom count {} != last frame {} in {}", n, *prev,
            filename);
        return IoStatus::InvalidArgument;
      }
    }
  }
 
  std::ofstream out;
  out.open(filename,
           append ? (std::ios::out | std::ios::app | std::ios::binary)
                  : (std::ios::out | std::ios::trunc | std::ios::binary));
  if (!out) {
    EONC_LOG_ERROR("matter2xyz: cannot open {}", filename);
    return IoStatus::OpenError;
  }
 
  // A frame header must start its own line. eOn's own frames end in a
  // newline; a hand-written or foreign tail may not.
  if (append && !ends_with_newline(filename)) {
    out.put('\n');
  }
 
  const Matrix3d cell = m.getCell();
  out << std::format(
      "{}\nLattice=\"{:.17g} {:.17g} {:.17g} {:.17g} {:.17g} {:.17g} "
      "{:.17g} {:.17g} {:.17g}\" Properties=species:S:1:pos:R:3 Generated "
      "by eOn\n",
      n, cell(0, 0), cell(0, 1), cell(0, 2), cell(1, 0), cell(1, 1), cell(1, 2),
      cell(2, 0), cell(2, 1), cell(2, 2));
  const AtomMatrix pos = m.getPositions();
  for (long i = 0; i < n; ++i) {
    out << std::format("{}\t{:.17g}\t{:.17g}\t{:.17g}\n",
                       symbol_for_z(m.getAtomicNr(i)), pos(i, 0), pos(i, 1),
                       pos(i, 2));
  }
  // Close before testing: the destructor's flush is where a full disk or a
  // short write surfaces, and by then the state is gone.
  out.close();
  if (!out) {
    EONC_LOG_ERROR("matter2xyz: failed to write {}", filename);
    return IoStatus::WriteError;
  }
  return IoStatus::Ok;
}
 
IoStatus writeTibble(Matter &m, std::string fname) {
  // Debug table, not a structure format. getForces()/getPotentialEnergy()
  // run computePotential() on a dirty pot, which would charge a dump to
  // the force-call count in results.dat. Cached values only; the header
  // drops the columns it cannot fill. atmID is the CON column-5 id.
  const bool have_forces = !m.needsForceUpdate();
  if (!have_forces) {
    EONC_LOG_WARNING(
        "writeTibble: pot is dirty, {} omits the force and energy columns",
        fname);
  }
  const AtomMatrix pos = m.getPositions();
  std::ofstream out(fname);
  if (!out) {
    EONC_LOG_ERROR("writeTibble: cannot open {}", fname);
    return IoStatus::OpenError;
  }
  out << (have_forces ? "x y z fx fy fz energy mass symbol atmID fixed\n"
                      : "x y z mass symbol atmID fixed\n");
  const AtomMatrix fSys = have_forces ? m.getForces() : AtomMatrix();
  const double eSys = have_forces ? m.getPotentialEnergy() : 0.0;
  for (long idx = 0; idx < m.numberOfAtoms(); ++idx) {
    out << std::format("{} {} {}", pos(idx, 0), pos(idx, 1), pos(idx, 2));
    if (have_forces) {
      out << std::format(" {} {} {} {}", fSys(idx, 0), fSys(idx, 1),
                         fSys(idx, 2), eSys);
    }
    const auto mask = m.getFixedMask(idx);
    const int fixed_bits =
        (mask[0] ? 1 : 0) | (mask[1] ? 2 : 0) | (mask[2] ? 4 : 0);
    out << std::format(" {} {} {} {}\n", m.getMass(idx),
                       symbol_for_z(m.getAtomicNr(idx)), m.getAtomIndex(idx),
                       fixed_bits);
  }
  out.close();
  if (!out) {
    EONC_LOG_ERROR("writeTibble: failed to write {}", fname);
    return IoStatus::WriteError;
  }
  return IoStatus::Ok;
}
 
std::vector<readcon::ConFrame>
buildNebPathFrames(const std::vector<std::shared_ptr<Matter>> &path,
                   const std::vector<ConFrameMetadata> &metadata_per_image) {
  std::vector<readcon::ConFrame> frames;
  if (path.empty() || path.size() != metadata_per_image.size()) {
    EONC_LOG_ERROR(
        "buildNebPathFrames: path/metadata size mismatch (path={}, meta={})",
        path.size(), metadata_per_image.size());
    return frames;
  }
  for (const auto &img : path) {
    if (!img) {
      EONC_LOG_ERROR("buildNebPathFrames: null Matter in path");
      return {};
    }
  }
 
  Matter &template_m = *path.front();
  template_m.applyPeriodicBoundaryIfEnabled();
 
  std::vector<uint64_t> atom_ids;
  std::array<std::string, 2> prebox;
  std::array<std::string, 2> postbox;
  collect_ids_headers(template_m, atom_ids, prebox, postbox);
 
  frames.reserve(path.size());
  try {
    auto seed = seed_builder(template_m, prebox, postbox, atom_ids);
    for (size_t i = 0; i < path.size(); ++i) {
      Matter &img = *path[i];
      img.applyPeriodicBoundaryIfEnabled();
      if (img.numberOfAtoms() != template_m.numberOfAtoms()) {
        EONC_LOG_ERROR(
            "buildNebPathFrames: image {} atom count {} != template {}", i,
            img.numberOfAtoms(), template_m.numberOfAtoms());
        return {};
      }
      auto builder = seed.clone();
      apply_frame_metadata(builder, &metadata_per_image[i]);
      apply_geometry(builder, img, /*with_velocities=*/false,
                     &metadata_per_image[i]);
      frames.push_back(builder.build());
    }
  } catch (const std::exception &e) {
    EONC_LOG_ERROR("buildNebPathFrames failed: {}", e.what());
    return {};
  }
  return frames;
}
 
IoStatus writeConFrames(std::string filename,
                        const std::vector<readcon::ConFrame> &frames) {
  if (frames.empty()) {
    return IoStatus::InvalidArgument;
  }
  filename = ensure_extension(std::move(filename), ".con");
  const fs::path path(filename);
  const std::lock_guard<std::mutex> guard(append_mutex());
  const auto key = append_key(path);
  const auto status = write_frames(path, frames, kConPrecision);
  if (io_ok(status)) {
    remember_stamp(key, path);
    mirror_con_corpus(path.string());
  } else {
    append_stamps().erase(key);
  }
  return status;
}
 
IoStatus writeNebPath(std::string filename,
                      const std::vector<std::shared_ptr<Matter>> &path,
                      const std::vector<ConFrameMetadata> &metadata_per_image) {
  auto frames = buildNebPathFrames(path, metadata_per_image);
  if (frames.empty()) {
    return IoStatus::InvalidArgument;
  }
  return writeConFrames(std::move(filename), frames);
}
 
} // namespace eonc::io