eonc::helpers Namespace Reference

Namespaces
namespace	create
namespace	eigen
namespace	neb_paths
namespace	surrogate

Functions
std::string	get_value_from_env_or_param (const char *env_variable, const std::string &param_value, const std::string &default_value, const std::string &warning_message, const bool is_mandatory)
std::unique_ptr< Job >	makeJob (std::unique_ptr< Parameters > params, Runtime &runtime)
	Borrow: caller keeps Runtime alive.
std::unique_ptr< Job >	makeJob (std::unique_ptr< Parameters > params, std::unique_ptr< Runtime > runtime)
	Own the unique_ptr Runtime.
std::unique_ptr< Job >	makeJob (std::unique_ptr< Parameters > params, Runtime &&runtime)
	Own the rvalue Runtime.
std::unique_ptr< Job >	makeJob (std::unique_ptr< Parameters > params)
	One-shot: Job owns a default Runtime.
std::shared_ptr< Potential >	makePotential (const Parameters &params)
std::shared_ptr< Potential >	makePotential (PotType ptype, const Parameters &params, Runtime &runtime)
std::shared_ptr< Potential >	makePotential (PotType ptype, const Parameters &params)
template<typename T>
std::vector< T >	get_val_from_string (std::string_view line, std::optional< size_t > nelements=std::nullopt)
	Parse a string into values.
template std::vector< size_t >	get_val_from_string (std::string_view, std::optional< size_t >)
std::vector< std::string >	get_split_strings (std::string_view line)
	Split a string into constituent strings.
bool	isNumber (std::string_view token)
	Figure out if a string has a number in it.
gpr::InputParameters	eon_parameters_to_gpr (const Parameters &parameters)
	Create a parameters object for gpr_dimer.
gpr::AtomsConfiguration	eon_matter_to_atmconf (Matter *matter)
	Create a configuration of atoms for gpr_dimer.
gpr::Observation	eon_matter_to_init_obs (Matter *matter)
	Create an initial Observation object for gpr_dimer.
AtomMatrix	makeOrthogonal (const AtomMatrix v1, const AtomMatrix v2)
bool	relaxMatter (Matter &matter, const Parameters &params, bool quiet=false, bool writeMovie=false, bool checkpoint=false, std::string prefixMovie=std::string(), std::string prefixCheckpoint=std::string(), std::vector< readcon::ConFrame > *outFrames=nullptr)
void	getTime (double *real, double *user, double *sys)
bool	existsFile (std::string filename)
std::optional< std::string >	enterJobDirectory (std::string_view jobPath)
	Enter jobPath as the working directory.
bool	stageReturnLog (std::string_view logHome, std::string_view name)
	Copy name from logHome into the current directory when that file is not already this directory's copy.
std::string	getRelevantFile (std::string filename)
VectorXd	loadMasses (std::string filename, int nAtoms)
AtomMatrix	loadMode (FILE *modeFile, int nAtoms)
AtomMatrix	loadMode (std::string filename, int nAtoms)
bool	loadOrSynthesizeDisplacement (Matter &target, const Matter &initial, const std::string &displacementPath, const std::string &modePath, double scale)
bool	applyClientDisplacement (Matter &target, const Matter &initial, const Parameters &params, AtomMatrix *modeOut)
void	saveMode (FILE *modeFile, std::shared_ptr< Matter > matter, AtomMatrix mode)
	Write a mode; constrained axes are emitted as 0.
void	saveMode (const std::string &filename, std::shared_ptr< Matter > matter, AtomMatrix mode)
std::vector< int >	split_string_int (std::string s, std::string delim)
std::optional< std::string_view >	convergenceMetricLabel (std::string_view metric)
	Display label for a force-convergence metric, or nullopt when the name is none of the four the optimizer-driven searches understand (norm, rms, max_atom, max_component).
void	requireKnownConvergenceMetric (std::string_view metric, std::string_view context)
	Throws std::invalid_argument naming context when metric is unrecognized.
template<typename T>
bool	eigenEquality (const Eigen::MatrixBase< T > &lhs, const Eigen::MatrixBase< T > &rhs, const double threshold=1e-4)
	Check two eigen objects for equality.

Variables
constexpr double	pi = 3.14159265358979323846

Function Documentation

applyClientDisplacement()

bool eonc::helpers::applyClientDisplacement	(	Matter &	target,
		const Matter &	initial,
		const Parameters &	params,
		AtomMatrix *	modeOut )

Definition at line 246 of file HelperFunctions.cpp.

                                                                 {
  using namespace eonc::EpiCenters;
  const auto &opt = params.saddle_search_options();
  const std::string &dtype = opt.displace_type;
  if (dtype == DISP_LOAD) {
    return false;
  }
 
  long epicenter = -1;
  const double cutoff = params.structure_comparison_options().neighbor_cutoff;
  if (dtype == DISP_LISTED_ATOMS) {
    epicenter = listedAtomEpiCenter(&initial, opt.displace_atom_list);
  } else if (dtype == DISP_RANDOM) {
    epicenter = randomFreeAtomEpiCenter(&initial);
  } else if (dtype == DISP_LAST_ATOM) {
    epicenter = lastAtom(&initial);
  } else if (dtype == DISP_MIN_COORDINATED) {
    epicenter = minCoordinatedEpiCenter(&initial, cutoff);
  } else if (dtype == DISP_NOT_FCC_OR_HCP) {
    epicenter = cnaEpiCenter(&initial, cutoff);
  } else {
    return false;
  }
 
  target = initial;
  const long n = initial.numberOfAtoms();
  const double radius = opt.displace_radius;
  const double mag = opt.displace_magnitude;
  AtomMatrix pos = initial.getPositionsCopy();
  AtomMatrix mode = AtomMatrix::Zero(n, 3);
  for (long i = 0; i < n; ++i) {
    if (initial.getFixed(i)) {
      continue;
    }
    const double dist = (i == epicenter) ? 0.0 : initial.distance(epicenter, i);
    if (dist <= radius) {
      for (int a = 0; a < 3; ++a) {
        mode(i, a) = eonc::rng::gaussRandom(0.0, mag);
      }
    }
  }
  const double norm = mode.norm();
  if (norm > 0.0) {
    pos += mode;
    mode /= norm;
  } else if (epicenter >= 0 && epicenter < n && !initial.getFixed(epicenter)) {
    mode(epicenter, 0) = 1.0;
    pos(epicenter, 0) += mag;
  }
  target.setPositions(pos);
  if (modeOut != nullptr) {
    *modeOut = std::move(mode);
  }
  return true;
}

convergenceMetricLabel()

std::optional< std::string_view > eonc::helpers::convergenceMetricLabel

(

std::string_view

metric

)

Display label for a force-convergence metric, or nullopt when the name is none of the four the optimizer-driven searches understand (norm, rms, max_atom, max_component).

Single source of the accepted spellings.

Definition at line 355 of file HelperFunctions.cpp.

                                                         {
  if (metric == "max_atom") {
    return "Max atom force";
  }
  if (metric == "max_component") {
    return "Max force comp";
  }
  if (metric == "norm") {
    return "||Force||";
  }
  if (metric == "rms") {
    return "RMS force";
  }
  return std::nullopt;
}

eigenEquality()

template<typename T>

bool eonc::helpers::eigenEquality	(	const Eigen::MatrixBase< T > &	lhs,
		const Eigen::MatrixBase< T > &	rhs,
		const double	threshold = 1e-4 )

Check two eigen objects for equality.

Parameters

Eigen::MatrixBase

the underlying base class

Definition at line 26 of file MatrixHelpers.hpp.

                                                {
  return lhs.isApprox(rhs, threshold);
}

enterJobDirectory()

std::optional< std::string > eonc::helpers::enterJobDirectory

(

std::string_view

jobPath

)

Enter jobPath as the working directory.

Returns nullopt on success. On failure the working directory is unchanged and the string is the filesystem error, so the caller must not run the job or report the path as finished.

Definition at line 80 of file HelperFunctions.cpp.

                                                     {
  try {
    std::filesystem::current_path(std::filesystem::path{jobPath});
  } catch (const std::filesystem::filesystem_error &err) {
    return std::string{err.what()};
  }
  return std::nullopt;
}

eon_matter_to_atmconf()

gpr::AtomsConfiguration eonc::helpers::eon_matter_to_atmconf

(

Matter *

matter

)

Create a configuration of atoms for gpr_dimer.

Parameters

*Matter

An eOn Matter object

Remember that the atom type in eOn is the real atomic number, while in GPR Dimer it is a set of values from 0 to n-1 so this is eOn

False "atomtype" for GPR Dimer

Is a fixed atom Use eOn's atomtype as a key for the GPR's fake atomtype

Is moving

Special case when there's only one atom type, we can now just use the set function of the Field. Essentially now we only have one atom type

End case where we have both nonzero moving and nonzero frozen atoms

Now we will consider the case when everything is moving Everything is almost exactly the same, only we don't have frozen atoms

Special case when there's only one atom type, we can now just use the set function of the Field. Essentially now we only have one atom type

Definition at line 112 of file GPRHelpers.cpp.

                                                                       {
  if (!matter) {
    throw std::invalid_argument("eon_matter_to_atmconf: null Matter");
  }
  gpr::AtomsConfiguration atoms_config;
  aux::ProblemSetUp problem_setup;
  gpr::Index_t number_of_mov_atoms;
  gpr::Index_t number_of_fro_atoms;
  std::set<int> unique_atomtypes;
  gpr::Index_t n_at;
  std::vector<int> atomnrs;
  std::unordered_map<int, int>
      atype_to_gprd_atype; 
  int fake_atype;          
 
  atoms_config.clear();
  // gpr_optim stores positions as row-major 1×(3N), matching AtomMatrix flat
  // layout.
  copyAtomMatrixToCoord(matter->getPositions(), atoms_config.positions);
  const auto nAtoms = matter->numberOfAtoms();
  atoms_config.is_frozen.resize(1, nAtoms);
  atoms_config.id.resize(1, nAtoms);
  atoms_config.atomicNrs.resize(1, nAtoms);
  for (auto i = 0; i < nAtoms; i++) {
    atomnrs.push_back(matter->getAtomicNr(i));
    // Field matrices are 1×N; use (0, i). getFixed is bool → MOVING/FROZEN.
    atoms_config.atomicNrs(0, i) = matter->getAtomicNr(i);
    atoms_config.is_frozen(0, i) =
        matter->getFixed(i) ? FROZEN_ATOM : MOVING_ATOM;
    atoms_config.id(0, i) = static_cast<gpr::Index_t>(i + 1);
  }
 
  unique_atomtypes = std::set<int>(atomnrs.begin(), atomnrs.end());
  n_at = unique_atomtypes.size();
  fake_atype = 0;
  for (auto uatom : unique_atomtypes) {
    atype_to_gprd_atype.insert(
        std::pair<int, int>(static_cast<int>(uatom), fake_atype));
    fake_atype++;
  }
 
  number_of_mov_atoms = atoms_config.countMovingAtoms();
  number_of_fro_atoms =
      static_cast<gpr::Index_t>(atoms_config.is_frozen.size()) -
      number_of_mov_atoms;
 
  if (number_of_fro_atoms > 0 && number_of_mov_atoms > 0) {
    // Resize structures for moving and frozen atoms
    atoms_config.atoms_mov.resize(number_of_mov_atoms);
    atoms_config.atoms_froz_inactive.resize(number_of_fro_atoms);
 
    // Dense GPR types from the Z→0..n-1 map.
    if (atype_to_gprd_atype.size() > 1) {
      int mov_counter = 0;
      int froz_inactive_counter = 0;
      for (auto i = 0; i < nAtoms; i++) {
        if (matter->getFixed(i)) {
          atoms_config.atoms_froz_inactive.type(0, froz_inactive_counter) =
              atype_to_gprd_atype.at(atomnrs[i]);
          froz_inactive_counter++;
        } else {
          atoms_config.atoms_mov.type(0, mov_counter) =
              atype_to_gprd_atype.at(atomnrs[i]);
          mov_counter++;
        }
      }
    }
    else if (atype_to_gprd_atype.size() == 1) {
      atoms_config.atoms_mov.type.setConstant(
          atype_to_gprd_atype.at(atomnrs[0]));
      atoms_config.atoms_froz_inactive.type.setConstant(
          atype_to_gprd_atype.at(atomnrs[0]));
    }
    // Assign moving and frozen atoms and list all frozen atoms as inactive
    gpr::Index_t counter_f = 0, counter_m = 0;
    for (gpr::Index_t n = 0;
         n < static_cast<gpr::Index_t>(atoms_config.is_frozen.size()); ++n) {
      if (atoms_config.is_frozen(0, n) == MOVING_ATOM)
        gpr::coord::set(atoms_config.atoms_mov.positions, 0, counter_m++,
                        gpr::coord::at(atoms_config.positions, n));
      else
        gpr::coord::set(atoms_config.atoms_froz_inactive.positions, 0,
                        counter_f++, gpr::coord::at(atoms_config.positions, n));
    }
  } else {
    if (number_of_mov_atoms == 0) {
      throw std::runtime_error("eon_matter_to_atmconf: no moving atoms");
    }
    atoms_config.atoms_mov.resize(number_of_mov_atoms);
 
    // Dense GPR types from the Z→0..n-1 map.
    if (atype_to_gprd_atype.size() > 1) {
      int mov_counter = 0;
      for (auto i = 0; i < nAtoms; i++) {
        atoms_config.atoms_mov.type(0, mov_counter) =
            atype_to_gprd_atype.at(atomnrs[i]);
        mov_counter++;
      }
    }
    else if (atype_to_gprd_atype.size() == 1) {
      atoms_config.atoms_mov.type.setConstant(
          atype_to_gprd_atype.at(atomnrs[0]));
    }
  }
  // Pairtype indices for pairs of atomtypes (n_at x n_at)
  // Active pairtypes are indexed as 0,1,...,n_pt-1. Inactive pairtypes are
  // given index EMPTY.
  atoms_config.pairtype.resize(n_at, n_at);
  atoms_config.pairtype.setConstant(EMPTY);
 
  // Set pairtype indices for moving+moving atom pairs (and update number of
  // active pairtypes)
  problem_setup.setPairtypeForMovingAtoms(
      atoms_config.atoms_mov.type, atoms_config.n_pt, atoms_config.pairtype);
 
  return atoms_config;
}

eon_matter_to_init_obs()

gpr::Observation eonc::helpers::eon_matter_to_init_obs

(

Matter *

matter

)

Create an initial Observation object for gpr_dimer.

Note that this is essentially only for the setup, and it does not actually append or handle the Observation structure other than for initialization of atomic gp dimer

Parameters

*Matter

An eOn Matter object

eon_parameters_to_gpr()

gpr::InputParameters eonc::helpers::eon_parameters_to_gpr

(

const Parameters &

parameters

)

Create a parameters object for gpr_dimer.

TODO: Get the cell size from a Matter object

Parameters

*parameters

An eOn parameters object

Definition at line 26 of file GPRHelpers.cpp.

                                                               {
  gpr::InputParameters p;
  // Problem parameters
  p.actdist_fro.value = parameters.gpr_dimer_options().active_radius;
  p.dimer_sep.value = parameters.gpr_dimer_options().dimer_sep;
  p.method_rot.value = parameters.gpr_dimer_options().rot_opt_method;
  p.method_trans.value = parameters.gpr_dimer_options().trans_opt_method;
  p.param_trans.value[0] = parameters.gpr_dimer_options().conv_step;
  p.param_trans.value[1] = parameters.gpr_dimer_options().max_step;
  // Saddle point convergence parameter
  p.T_dimer.value = parameters.saddle_search_options().converged_force;
  p.T_anglerot_init.value = parameters.gpr_dimer_options().converged_angle;
  // ---
  p.initrot_nogp.value = parameters.gpr_dimer_options().init_rot_gp;
  p.num_iter_initrot.value = parameters.gpr_dimer_options().init_rotations_max;
  // unused except in the capnp for now
  // p.inittrans_nogp.value = parameters.gpr_dimer_options().init_trans_gp;
  p.T_anglerot_gp.value = parameters.gpr_dimer_options().relax_conv_angle;
  p.num_iter_rot_gp.value = parameters.gpr_dimer_options().relax_rotations_max;
  p.divisor_T_dimer_gp.value =
      parameters.gpr_dimer_options().divisor_t_dimer_gp;
  p.disp_max.value = parameters.gpr_dimer_options().midpoint_max_disp;
  p.ratio_at_limit.value = parameters.gpr_dimer_options().ratio_at_limit;
  p.num_bigiter.value = parameters.gpr_dimer_options().max_outer_iterations;
  p.num_iter.value = parameters.gpr_dimer_options().max_inner_iterations;
  p.islarge_num_iter.value = parameters.gpr_dimer_options().many_iterations;
  // GPR Parameters
  p.gp_sigma2.value = parameters.gpr_dimer_options().gpr_params.sigma2;
  p.jitter_sigma2.value =
      parameters.gpr_dimer_options().gpr_params.jitter_sigma2;
  p.sigma2.value = parameters.gpr_dimer_options().gpr_params.noise_sigma2;
  p.prior_mu.value = parameters.gpr_dimer_options().gpr_params.prior_mu;
  p.prior_nu.value = parameters.gpr_dimer_options().gpr_params.prior_nu;
  p.prior_s2.value = parameters.gpr_dimer_options().gpr_params.prior_sigma2;
  // gpr_optim enables the check only when this string equals "true".
  // A bool assigns as one char, so the flag never turns on.
  p.check_derivative.value =
      parameters.gpr_dimer_options().opt_params.check_derivatives ? "true"
                                                                  : "false";
  p.max_iter.value = parameters.gpr_dimer_options().opt_params.max_iterations;
  p.tolerance_func.value = parameters.gpr_dimer_options().opt_params.tol_func;
  p.tolerance_sol.value = parameters.gpr_dimer_options().opt_params.tol_sol;
  p.lambda_limit.value = parameters.gpr_dimer_options().opt_params.lambda_limit;
  p.lambda.value = parameters.gpr_dimer_options().opt_params.lambda_init;
  // Prune
  p.use_prune.value = parameters.gpr_dimer_options().prune_params.use_prune;
  p.start_prune_at.value = parameters.gpr_dimer_options().prune_params.begin;
  p.nprune_vals.value = parameters.gpr_dimer_options().prune_params.n_vals;
  p.prune_threshold.value =
      parameters.gpr_dimer_options().prune_params.threshold;
  // Debugging
  p.report_level.value =
      parameters.gpr_dimer_options().debug_params.report_level;
  p.debug_level.value = parameters.gpr_dimer_options().debug_params.debug_level;
  p.debug_output_dir.value =
      parameters.gpr_dimer_options().debug_params.out_dir;
  p.debug_output_file_R.value =
      parameters.gpr_dimer_options().debug_params.pos_file;
  p.debug_output_file_E.value =
      parameters.gpr_dimer_options().debug_params.energy_file;
  p.debug_output_file_G.value =
      parameters.gpr_dimer_options().debug_params.grad_file;
  p.debug_output_file_extension.value =
      parameters.gpr_dimer_options().debug_params.out_ext;
  p.debug_offset_from_mid_point.value =
      parameters.gpr_dimer_options().debug_params.offset_mid_point;
  p.debug_dy.value = parameters.gpr_dimer_options().debug_params.dy;
  p.debug_dz.value = parameters.gpr_dimer_options().debug_params.dz;
  return p;
}

existsFile()

bool eonc::helpers::existsFile

(

std::string

filename

)

Definition at line 75 of file HelperFunctions.cpp.

                                               {
  return std::filesystem::exists(filename);
}

get_split_strings()

std::vector< std::string > eonc::helpers::get_split_strings

(

std::string_view

line

)

Split a string into constituent strings.

Based on https://www.fluentcpp.com/2017/04/21/how-to-split-a-string-in-c/

Parameters

line	A thing to be parsed

Definition at line 46 of file StringHelpers.cc.

                                                            {
  std::istringstream ss{std::string{line}};
  std::vector<std::string> split_strings{std::istream_iterator<std::string>{ss},
                                         std::istream_iterator<std::string>()};
  return split_strings;
}

get_val_from_string() [1/2]

template<typename T>

std::vector< T > eonc::helpers::get_val_from_string	(	std::string_view	line,
		std::optional< size_t >	nelements = std::nullopt )

Parse a string into values.

Parameters

line	A thing to be parsed

Definition at line 9 of file StringHelpers.cc.

                                                                  {
  assert(not line.empty());
  std::vector<T> retval;
  const bool b_isunsigned{std::is_unsigned<T>::value};
  auto elements{get_split_strings(line)};
  if (nelements.has_value()) {
    // Used to truncate if the number of elements is given
    assert(nelements > 0);
    elements.resize(nelements.value());
  }
  // If it is unsigned then use long double else T
  for (typename std::conditional<b_isunsigned, long double, T>::type tmp;
       auto elem : elements) {
    if (not isNumber(elem)) {
      continue;
    }
    std::istringstream ss{
        elem}; // instead of {ss.str(elem); ss >> tmp; ss.clear();}
    ss >> tmp;
    if (b_isunsigned and tmp < 0) {
      std::cerr << "Can't represent negative numbers with an unsigned type, "
                   "bailing on "
                << tmp << "\n";
      assert(tmp > 0);
    }
    retval.push_back(tmp);
  }
  return retval;
}

get_val_from_string() [2/2]

template std::vector< size_t > eonc::helpers::get_val_from_string	(	std::string_view	,
		std::optional< size_t >	)

get_value_from_env_or_param()

std::string eonc::helpers::get_value_from_env_or_param ( const char * env_variable, const std::string & param_value, const std::string & default_value, const std::string & warning_message, const bool is_mandatory )

nodiscard

Definition at line 7 of file EnvHelpers.cc.

                                                                 {
  const char *env_value = std::getenv(env_variable);
  if (env_value != nullptr) {
    return std::string(env_value);
  }
 
  if (!param_value.empty()) {
    return param_value;
  }
 
  if (is_mandatory) {
    throw std::runtime_error(
        "Environment variable " + std::string(env_variable) +
        " is not set and no parameter value provided. Please set it in the "
        "configuration or as an environment variable.\n");
  }
 
  if (!default_value.empty() && !warning_message.empty()) {
    EONC_LOG_WARNING("{}", warning_message);
  }
 
  return default_value;
}

getRelevantFile()

std::string eonc::helpers::getRelevantFile

(

std::string

filename

)

Definition at line 122 of file HelperFunctions.cpp.

                                                         {
  const auto dot = filename.rfind('.');
  const std::string prefix =
      (dot == std::string::npos) ? filename : filename.substr(0, dot);
  const std::string postfix =
      (dot == std::string::npos) ? std::string{} : filename.substr(dot);
  std::string filenameRelevant = prefix + "_cp" + postfix;
  if (existsFile(filenameRelevant)) {
    return filenameRelevant;
  }
  filenameRelevant = prefix + "_in" + postfix;
  if (existsFile(filenameRelevant)) {
    return filenameRelevant;
  }
  return filename;
}

getTime()

void eonc::helpers::getTime	(	double *	real,
		double *	user,
		double *	sys )

Definition at line 47 of file HelperFunctions.cpp.

                                                                 {
  using namespace std::chrono;
  auto now = steady_clock::now();
  if (real) {
    *real = duration<double>(now.time_since_epoch()).count();
  }
 
#ifdef _WIN32
  if (user)
    *user = 0.0;
  if (sys)
    *sys = 0.0;
#else
  struct rusage r_usage;
  if (getrusage(RUSAGE_SELF, &r_usage) != 0) {
    EONC_LOG_WARNING("problem getting usage info: {}", strerror(errno));
  }
  if (user) {
    *user = static_cast<double>(r_usage.ru_utime.tv_sec) +
            static_cast<double>(r_usage.ru_utime.tv_usec) / 1e6;
  }
  if (sys) {
    *sys = static_cast<double>(r_usage.ru_stime.tv_sec) +
           static_cast<double>(r_usage.ru_stime.tv_usec) / 1e6;
  }
#endif
}

isNumber()

bool eonc::helpers::isNumber

(

std::string_view

token

)

Figure out if a string has a number in it.

Based on

Definition at line 53 of file StringHelpers.cc.

                                    {
  return std::regex_match(
      std::string{token},
      std::regex(("((\\+|-)?[[:digit:]]+)(\\.(([[:digit:]]+)?))?")));
}

loadMasses()

VectorXd eonc::helpers::loadMasses	(	std::string	filename,
		int	nAtoms )

Definition at line 139 of file HelperFunctions.cpp.

                                                               {
  std::ifstream massFile(filename.c_str());
  if (!massFile.is_open()) {
    EONC_LOG_CRITICAL("File {} was not found", filename);
    throw std::runtime_error(std::format("cannot open {}", filename));
  }
 
  VectorXd masses(nAtoms);
  for (int i = 0; i < nAtoms; i++) {
    double mass;
    if (!(massFile >> mass)) {
      EONC_LOG_CRITICAL("Error reading {}", filename);
      throw std::runtime_error(
          std::format("{} ended after {} of {} masses", filename, i, nAtoms));
    }
    masses(i) = mass;
  }
 
  massFile.close();
 
  return masses;
}

loadMode() [1/2]

AtomMatrix eonc::helpers::loadMode	(	FILE *	modeFile,
		int	nAtoms )

Definition at line 162 of file HelperFunctions.cpp.

                                                           {
  AtomMatrix mode;
  mode.resize(nAtoms, 3);
  mode.setZero();
  for (int i = 0; i < nAtoms; i++) {
    if (fscanf(modeFile, "%lf %lf %lf", &mode(i, 0), &mode(i, 1),
               &mode(i, 2)) != 3) {
      EONC_LOG_CRITICAL("Mode file ended after {} of {} atoms", i, nAtoms);
      throw std::runtime_error(
          std::format("mode file ended after {} of {} atoms", i, nAtoms));
    }
  }
  return mode;
}

loadMode() [2/2]

AtomMatrix eonc::helpers::loadMode	(	std::string	filename,
		int	nAtoms )

Definition at line 177 of file HelperFunctions.cpp.

                                                               {
  // Unique FILE* with RAII cleanup
  auto closer = [](FILE *f) {
    if (f)
      std::fclose(f);
  };
  std::unique_ptr<FILE, decltype(closer)> modeFile(
      std::fopen(filename.c_str(), "rb"), closer);
  if (!modeFile) {
    EONC_LOG_CRITICAL("File {} was not found", filename);
    throw std::runtime_error(std::format("cannot open {}", filename));
  }
  return loadMode(modeFile.get(), nAtoms);
}

loadOrSynthesizeDisplacement()

bool eonc::helpers::loadOrSynthesizeDisplacement	(	Matter &	target,
		const Matter &	initial,
		const std::string &	displacementPath,
		const std::string &	modePath,
		double	scale )

Definition at line 192 of file HelperFunctions.cpp.

                                             {
  if (eonc::io::io_ok(target.con2matter(displacementPath))) {
    if (target.numberOfAtoms() != initial.numberOfAtoms()) {
      EONC_LOG_ERROR("{} holds {} atoms, the initial structure has {}",
                     displacementPath, target.numberOfAtoms(),
                     initial.numberOfAtoms());
      return false;
    }
    // displacement.con may carry stale fixed-atom coordinates from a prior run.
    // It also usually has sequential column-5 ids; keep the reactant's.
    const AtomMatrix &initPos = initial.getPositions();
    AtomMatrix pos = target.getPositionsCopy();
    const long n = initial.numberOfAtoms();
    std::vector<long> fileMap(static_cast<size_t>(n));
    for (long i = 0; i < n; i++) {
      if (initial.getFixed(i)) {
        pos.row(i) = initPos.row(i);
      }
      target.setAtomIndex(i, initial.getAtomIndex(i));
      fileMap[static_cast<size_t>(i)] = initial.mapFileRow(i);
    }
    target.setFileToMatter(std::move(fileMap));
    target.setPositions(pos);
    return true;
  }
  if (!existsFile(modePath)) {
    return false;
  }
  AtomMatrix mode =
      loadMode(modePath, static_cast<int>(initial.numberOfAtoms()));
  const double norm = mode.norm();
  if (!(norm > 0.0)) {
    return false;
  }
  mode *= (scale / norm);
  target = initial;
  AtomMatrix pos = initial.getPositionsCopy();
  pos += mode;
  const AtomMatrix &initPos = initial.getPositions();
  const long n = initial.numberOfAtoms();
  for (long i = 0; i < n; i++) {
    if (initial.getFixed(i)) {
      pos.row(i) = initPos.row(i);
    }
  }
  target.setPositions(pos);
  EONC_LOG_INFO("Synthesized displacement from pos.con + scale {:.6g} * unit "
                "mode in {} (missing {})",
                scale, modePath, displacementPath);
  return true;
}

makeJob() [1/4]

std::unique_ptr< Job > eonc::helpers::makeJob

(

std::unique_ptr< Parameters >

params

)

One-shot: Job owns a default Runtime.

Definition at line 65 of file JobRegistry.cpp.

                                                             {
  return makeJob(std::move(params), std::make_unique<Runtime>());
}

makeJob() [2/4]

std::unique_ptr< Job > eonc::helpers::makeJob	(	std::unique_ptr< Parameters >	params,
		Runtime &&	runtime )

Own the rvalue Runtime.

Definition at line 59 of file JobRegistry.cpp.

                                                {
  return makeJob(std::move(params),
                 std::make_unique<Runtime>(std::move(runtime)));
}

makeJob() [3/4]

std::unique_ptr< Job > eonc::helpers::makeJob	(	std::unique_ptr< Parameters >	params,
		Runtime &	runtime )

Borrow: caller keeps Runtime alive.

Definition at line 43 of file JobRegistry.cpp.

                                               {
  return makeJobFromFactory(std::move(params), runtime);
}

makeJob() [4/4]

std::unique_ptr< Job > eonc::helpers::makeJob	(	std::unique_ptr< Parameters >	params,
		std::unique_ptr< Runtime >	runtime )

Own the unique_ptr Runtime.

Definition at line 48 of file JobRegistry.cpp.

                                                             {
  if (!runtime) {
    throw std::logic_error("makeJob: null Runtime");
  }
  Runtime &ref = *runtime;
  auto job = makeJobFromFactory(std::move(params), ref);
  job->adoptRuntime(std::move(runtime));
  return job;
}

makeOrthogonal()

AtomMatrix eonc::helpers::makeOrthogonal	(	const AtomMatrix	v1,
		const AtomMatrix	v2 )

Definition at line 42 of file HelperFunctions.cpp.

                                                              {
  return v1 - matDot(v1, v2) * eonc::safemath::safe_normalized(v2);
}

makePotential() [1/3]

std::shared_ptr< Potential > eonc::helpers::makePotential

(

const Parameters &

params

)

Definition at line 246 of file Potential.cpp.

                                                                 {
  return makePotential(params.potential_options().potential, params);
}

makePotential() [2/3]

std::shared_ptr< Potential > eonc::helpers::makePotential	(	PotType	ptype,
		const Parameters &	params )

Definition at line 255 of file Potential.cpp.

                                                                   {
  // Inject config-file path before any potential constructor runs.
  // Called on every code path including Job::Job which uses this overload.
  PluginLoader::instance().add_config_paths(
      params.potential_options().potentialsPath);
  switch (ptype) {
  // TODO: Every potential must know their own type
  case PotType::EMT: {
    return (std::make_shared<EffectiveMediumTheory>(params));
    break;
  }
  case PotType::EXT_POT: {
    return (std::make_shared<ExtPot>(params));
    break;
  }
  case PotType::LJ: {
    return makeRgpot<rgpot::LJPot>(PotType::LJ, params, rgpot::LJConfig{});
    break;
  }
  case PotType::LJCLUSTER: {
    return makeRgpot<rgpot::LJClusterPot>(PotType::LJCLUSTER, params,
                                          rgpot::LJClusterConfig{});
    break;
  }
  case PotType::MORSE_PT: {
    return makeRgpot<rgpot::MorsePot>(PotType::MORSE_PT, params,
                                      rgpot::MorseConfig{});
    break;
  }
#ifdef CUH2_POT
  case PotType::CUH2: {
    return makeRgpotDefault<rgpot::fortranpots::CuH2Pot>(PotType::CUH2, params);
    break;
  }
#endif
#ifdef WITH_WATER
  case PotType::TIP4P: {
    return (std::make_shared<Tip4p>(params));
    break;
  }
  case PotType::SPCE: {
    return (std::make_shared<SpceCcl>(params));
    break;
  }
#ifdef WITH_FORTRAN
  case PotType::TIP4P_PT: {
    return (std::make_shared<Tip4p_Pt>(params));
    break;
  }
  case PotType::TIP4P_H: {
    return makeRgpotDefault<rgpot::fortranpots::WaterHPot>(PotType::TIP4P_H,
                                                           params);
    break;
  }
#endif
#endif
  // Fortran potentials: always available, loaded at runtime via dlopen
  case PotType::EAM_AL: {
    return makeRgpotDefault<rgpot::fortranpots::EAMAlPot>(PotType::EAM_AL,
                                                          params);
    break;
  }
  case PotType::EDIP: {
    return makeRgpotDefault<rgpot::fortranpots::EDIPPot>(PotType::EDIP, params);
    break;
  }
  case PotType::FEHE: {
    return makeRgpotDefault<rgpot::fortranpots::FeHePot>(PotType::FEHE, params);
    break;
  }
  case PotType::LENOSKY_SI: {
    return makeRgpotDefault<rgpot::fortranpots::LenoskyPot>(PotType::LENOSKY_SI,
                                                            params);
    break;
  }
  case PotType::SW_SI: {
    return makeRgpotDefault<rgpot::fortranpots::SWPot>(PotType::SW_SI, params);
    break;
  }
  case PotType::TERSOFF_SI: {
    return makeRgpotDefault<rgpot::fortranpots::TersoffPot>(PotType::TERSOFF_SI,
                                                            params);
    break;
  }
#ifndef _WIN32
#ifdef WITH_VASP
  case PotType::VASP: {
    return (std::make_shared<VASP>(params));
    break;
  }
#endif
#endif
  case PotType::LAMMPS: {
    return std::make_shared<LAMMPSPot>(params);
  }
#ifdef EONMPI
  case PotType::MPI: {
    return (std::make_shared<MPIPot>(params));
    break;
  }
#endif
#ifdef EMBED_PYTHON
#ifdef WITH_ASE_POT
  case PotType::ASE_POT: {
    return (std::make_shared<ASE>(params));
    break;
  }
#endif
#endif
#ifdef WITH_AMS
  case PotType::AMS: {
    return (std::make_shared<AMS>(params));
    break;
  }
  case PotType::AMS_IO: {
    return (std::make_shared<AMS_IO>(params));
    break;
  }
#endif
#ifdef WITH_CATLEARN
  case PotType::CatLearn: {
    return (std::make_shared<CatLearnPot>(params));
    break;
  }
#endif
// TODO: Handle Fortran interaction
#ifdef WITH_XTB
  case PotType::XTB: {
    return (std::make_shared<XTBPot>(params));
    break;
  }
#endif
#ifdef WITH_ASE_ORCA
  case PotType::ASE_ORCA: {
    return (std::make_shared<ASEOrcaPot>(params));
    break;
  }
#endif
#ifdef WITH_ASE_NWCHEM
  case PotType::ASE_NWCHEM: {
    return (std::make_shared<ASENwchemPot>(params));
    break;
  }
#endif
#ifdef WITH_METATOMIC
  case PotType::METATOMIC: {
    return (std::make_shared<MetatomicPotential>(params));
    break;
  }
#endif
#ifdef RGPOT_HAS_DFTD3
  case PotType::DFTD3: {
    rgpot::D3Damping damp = rgpot::D3Damping::BJ;
    std::string dname = params.dftd_options().d3_damping;
    for (char &c : dname) {
      c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
    }
    if (dname == "zero") {
      damp = rgpot::D3Damping::Zero;
    }
    return makeRgpot<rgpot::D3Pot>(
        PotType::DFTD3, params,
        rgpot::D3Config{.damping = damp,
                        .functional = params.dftd_options().functional,
                        .atm = params.dftd_options().atm});
    break;
  }
#endif
#ifdef RGPOT_HAS_DFTD4
  case PotType::DFTD4: {
    return makeRgpot<rgpot::D4Pot>(
        PotType::DFTD4, params,
        rgpot::D4Config{.functional = params.dftd_options().functional,
                        .charge = params.dftd_options().d4_charge,
                        .atm = params.dftd_options().atm});
    break;
  }
#endif
  case PotType::ZBL: {
    return makeRgpot<rgpot::ZBLPot>(
        PotType::ZBL, params,
        rgpot::ZBLConfig{
            .cut_inner = params.zbl_options().cut_inner,
            .cut_global = params.zbl_options().cut_global,
        });
    break;
  }
#ifndef IS_WINDOWS
  case PotType::SocketNWChem: {
    return (std::make_shared<SocketNWChemPot>(params));
    break;
  }
#endif
#ifdef WITH_RGPOT
  case PotType::RGPOT: {
    return (std::make_shared<RgpotPot>(params));
    break;
  }
#endif
  case PotType::MOPAC: {
    return makeRgpot<rgpot::MOPACPot>(
        PotType::MOPAC, params,
        rgpot::MOPACPot::Config{
            .charge = params.mopac_options().charge,
            .spin = params.mopac_options().spin,
            .model = params.mopac_options().model,
            .engine_path = params.mopac_options().engine_path,
        });
    break;
  }
#ifdef RGPOT_HAS_EXPR
  case PotType::EXPR: {
    if (params.expr_options().expression.empty()) {
      throw std::runtime_error(
          "ExprPot needs [ExprPot] expression, e.g. 0.5*lj + d3");
    }
    rgpot::ExprPot expr(params.expr_options().expression,
                        parse_expr_terms(params));
    return std::make_shared<RgpotAdapter<rgpot::ExprPot>>(PotType::EXPR, params,
                                                          std::move(expr));
    break;
  }
#endif
  default:
    EONC_LOG_ERROR("No known potential could be constructed from {}",
                   magic_enum::enum_name(ptype));
    eonc::log::get()->flush_log();
    throw std::runtime_error("Terminating");
    break;
  }
}

makePotential() [3/3]

std::shared_ptr< Potential > eonc::helpers::makePotential	(	PotType	ptype,
		const Parameters &	params,
		Runtime &	runtime )

Definition at line 250 of file Potential.cpp.

                                                                         {
  PotentialConstructionScope scope(runtime.pots());
  runtime.plugins().add_config_paths(params.potential_options().potentialsPath);
  return makePotential(ptype, params);
}

relaxMatter()

bool eonc::helpers::relaxMatter	(	Matter &	matter,
		const Parameters &	params,
		bool	quiet = false,
		bool	writeMovie = false,
		bool	checkpoint = false,
		std::string	prefixMovie = std::string(),
		std::string	prefixCheckpoint = std::string(),
		std::vector< readcon::ConFrame > *	outFrames = nullptr )

Definition at line 447 of file HelperFunctions.cpp.

                                                                       {
  eonc::log::Scoped m_log;
  auto objf = std::make_shared<MatterObjectiveFunction>(matter, params);
  auto optim = eonc::helpers::create::mkOptim(
      objf, params.optimizer_options().method, params);
 
  std::ostringstream min;
  min << prefixMovie;
  std::string minDatFilename = prefixMovie + ".dat";
  auto write_movie_frame = [&](uint64_t frameIndex, bool append,
                               double stepSize) {
    eonc::io::ConFrameMetadata metadata;
    metadata.frame_index = frameIndex;
    metadata.energy = matter.getPotentialEnergy();
    metadata.scalars.push_back({"step_size", stepSize});
    metadata.scalars.push_back({"convergence", objf->getConvergence()});
    if (outFrames) {
      outFrames->push_back(eonc::io::matterToConFrame(matter, &metadata));
    }
    if (writeMovie) {
      if (!eonc::io::io_ok(matter.matter2con(min.str(), append, &metadata))) {
        QUILL_LOG_WARNING(m_log, "Failed to write movie frame {}", min.str());
      }
    }
 
    if (params.debug_options().write_deprecated_outs) {
      std::ofstream minDat(minDatFilename,
                           append ? (std::ios::binary | std::ios::app)
                                  : std::ios::binary);
      if (minDat) {
        if (!append) {
          minDat << "iteration\tstep_size\tconvergence\tenergy\n";
        }
        minDat << std::format("{}\t{:.5e}\t{:.5e}\t{:.6f}\n", frameIndex,
                              stepSize, objf->getConvergence(),
                              matter.getPotentialEnergy());
      }
    }
  };
  if (writeMovie || outFrames) {
    write_movie_frame(0, false, 0.0);
  }
 
  int iteration = 0;
  if (!quiet) {
    QUILL_LOG_DEBUG(m_log, "{} {:10s}  {:14s}  {:18s}  {:13s}\n", "[Matter]",
                    "Iter", "Step size",
                    params.optimizer_options().convergence_metric_label,
                    "Energy");
    QUILL_LOG_DEBUG(m_log, "{} {:10}  {:14.5e}  {:18.5e}  {:13.5f}\n",
                    "[Matter]", iteration, 0.0, objf->getConvergence(),
                    matter.getPotentialEnergy());
  }
 
  while (!objf->isConverged() &&
         iteration < params.optimizer_options().max_iterations) {
 
    AtomMatrix pos = matter.getPositions();
 
    optim->step(params.optimizer_options().max_move);
    iteration++;
 
    double stepSize =
        eonc::geometry::maxAtomMotion(matter.pbc(matter.getPositions() - pos));
 
    if (!quiet) {
      QUILL_LOG_DEBUG(m_log, "{} {:10}  {:14.5e}  {:18.5e}  {:13.5f}",
                      "[Matter]", iteration, stepSize, objf->getConvergence(),
                      matter.getPotentialEnergy());
    }
 
    if (writeMovie || outFrames) {
      write_movie_frame(static_cast<uint64_t>(iteration), true, stepSize);
    }
 
    if (checkpoint) {
      std::ostringstream chk;
      chk << prefixCheckpoint << "_cp";
      if (!eonc::io::io_ok(matter.matter2con(chk.str(), false))) {
        QUILL_LOG_WARNING(m_log, "Failed to write checkpoint {}", chk.str());
      }
    }
  }
 
  if (iteration == 0) {
    if (!quiet) {
      QUILL_LOG_DEBUG(m_log, "{} {:10}  {:14.5e}  {:18.5e}  {:13.5f}",
                      "[Matter]", iteration, 0.0, objf->getConvergence(),
                      matter.getPotentialEnergy());
    }
  }
  return objf->isConverged();
}

requireKnownConvergenceMetric()

void eonc::helpers::requireKnownConvergenceMetric	(	std::string_view	metric,
		std::string_view	context )

Throws std::invalid_argument naming context when metric is unrecognized.

Call at construction: the metric is compared once per optimizer step, and a typo should not surface on iteration 4000.

Definition at line 371 of file HelperFunctions.cpp.

                                                                          {
  if (convergenceMetricLabel(metric)) {
    return;
  }
  throw std::invalid_argument(
      std::format("{} unknown convergence_metric: {}", context, metric));
}

saveMode() [1/2]

void eonc::helpers::saveMode	(	const std::string &	filename,
		std::shared_ptr< Matter >	matter,
		AtomMatrix	mode )

Definition at line 316 of file HelperFunctions.cpp.

                                                                            {
  std::ofstream out(filename);
  if (!out)
    return;
  const AtomMatrix free = matter->getFree();
  long const nAtoms = matter->numberOfAtoms();
  for (long i = 0; i < nAtoms; ++i) {
    out << std::format("{:.17g}\t{:.17g}\t{:.17g}\n", free(i, 0) * mode(i, 0),
                       free(i, 1) * mode(i, 1), free(i, 2) * mode(i, 2));
  }
}

saveMode() [2/2]

void eonc::helpers::saveMode	(	FILE *	modeFile,
		std::shared_ptr< Matter >	matter,
		AtomMatrix	mode )

Write a mode; constrained axes are emitted as 0.

Definition at line 305 of file HelperFunctions.cpp.

                                              {
  const AtomMatrix free = matter->getFree();
  long const nAtoms = matter->numberOfAtoms();
  for (long i = 0; i < nAtoms; ++i) {
    fprintf(modeFile, "%.17g\t%.17g\t%.17g\n", free(i, 0) * mode(i, 0),
            free(i, 1) * mode(i, 1), free(i, 2) * mode(i, 2));
  }
  return;
}

split_string_int()

std::vector< int > eonc::helpers::split_string_int	(	std::string	s,
		std::string	delim )

Definition at line 329 of file HelperFunctions.cpp.

                                                                  {
  std::vector<int> list;
  if (s.empty())
    return list;
 
  size_t start = 0;
  size_t end = s.find_first_of(delim);
  while (start < s.size()) {
    auto token = s.substr(start, end - start);
    if (!token.empty()) {
      try {
        list.push_back(std::stoi(token));
      } catch (const std::exception &) {
        return {}; // Parse error
      }
    }
    if (end == std::string::npos)
      break;
    start = end + 1;
    end = s.find_first_of(delim, start);
  }
  return list;
}

stageReturnLog()

bool eonc::helpers::stageReturnLog	(	std::string_view	logHome,
		std::string_view	name )

Copy name from logHome into the current directory when that file is not already this directory's copy.

True only when a regular file of that name is then present.

Definition at line 89 of file HelperFunctions.cpp.

                                                        {
  namespace fs = std::filesystem;
  if (name.empty() || name == "." || name == ".." ||
      name.find('/') != std::string_view::npos ||
      name.find('\\') != std::string_view::npos) {
    return false;
  }
  const fs::path dest{std::string{name}};
  std::error_code ec;
  const bool destFile = fs::is_regular_file(dest, ec);
  if (logHome.empty()) {
    return destFile;
  }
  const fs::path src = fs::path{std::string{logHome}} / dest.filename();
  ec.clear();
  if (!fs::is_regular_file(src, ec)) {
    return destFile;
  }
  ec.clear();
  if (destFile && fs::equivalent(src, dest, ec)) {
    return true;
  }
  ec.clear();
  fs::copy_file(src, dest, fs::copy_options::overwrite_existing, ec);
  if (ec) {
    EONC_LOG_ERROR("stageReturnLog: cannot copy {} to {}: {}", src.string(),
                   dest.string(), ec.message());
    return fs::is_regular_file(dest);
  }
  return true;
}

Variable Documentation

pi

double eonc::helpers::pi = 3.14159265358979323846

inlineconstexpr

Definition at line 29 of file HelperFunctions.h.