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)
std::shared_ptr< Potential >	makePotential (const Parameters &params)
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::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)
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 three the optimizer-driven searches understand.
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.
double	gaussRandom (double avg, double std)
double	random (long newSeed=0)
double	randomDouble ()
long	randomInt (int lower, int upper)
bool	identical (const Matter &m1, const Matter &m2, const double distanceDifference)
double	maxAtomMotion (const AtomMatrix v1)
AtomMatrix	maxAtomMotionApplied (const AtomMatrix v1, double maxMotion)
VectorXd	maxAtomMotionAppliedV (const VectorXd v1, double maxMotion)
double	maxAtomMotionV (const VectorXd v1)
AtomMatrix	maxMotionApplied (const AtomMatrix v1, double maxMotion)
VectorXd	maxMotionAppliedV (const VectorXd v1, double maxMotion)
long	numAtomsMoved (const AtomMatrix v1, double cutoff)
void	projectOutRotTrans (Eigen::VectorXd &step, const AtomMatrix &positions)
void	pushApart (std::shared_ptr< Matter > m1, double minDistance)
RotationMatrix	rotationExtract (const AtomMatrix r1, const AtomMatrix r2)
bool	rotationMatch (const Matter &m1, const Matter &m2, const double max_diff)
void	rotationRemove (const AtomMatrix r1, std::shared_ptr< Matter > m2)
bool	sortedR (const Matter &m1, const Matter &m2, const double distanceDifference)
void	translationRemove (Matter &m1, const AtomMatrix r1)

Variables
constexpr double	pi = 3.14159265358979323846

Function Documentation

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 three the optimizer-driven searches understand.

Single source of the accepted spellings, shared by the INI loader and every search that dispatches on the string.

Definition at line 250 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||";
  }
  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);
}

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

Does a horrible to ensure that this is filled correctly. Essentially we use the Map of <eOn atomtype, GPR faketype> to generate the fully filled vectors for moving and frozen_inactive FIXME: We should really just use the eOn atomtype everywhere

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

Sometimes, nothing happens

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

FIXME: Same caveats as documented above

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 100 of file GPRHelpers.cpp.

                                                                       {
  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);

    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) {
      QUILL_LOG_CRITICAL(
          eonc::log::get(),
          " You need to have atoms move!!!\nIn stillness there is only "
          "death\n");
      std::exit(1);
    }
    atoms_config.atoms_mov.resize(number_of_mov_atoms);

    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

Definition at line 235 of file GPRHelpers.cpp.

                                                                 {
  gpr::Observation o;
  o.clear();
  copyAtomMatrixToCoord(matter->getPositions(), o.R);
  // Forces as negative gradients; same N×3 → 1×(3N) packing as positions.
  AtomMatrix neg_forces = -matter->getForces();
  copyAtomMatrixToCoord(neg_forces, o.G);
  o.E.resize(1, 1);
  o.E(0, 0) = matter->getPotentialEnergy();
  return o;
}

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 23 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;
  p.check_derivative.value =
      parameters.gpr_dimer_options.opt_params.check_derivatives;
  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 72 of file HelperFunctions.cpp.

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

gaussRandom()

double eonc::rng::gaussRandom	(	double	avg,
		double	std )

Definition at line 41 of file RandomNumbers.cpp.

                                                  {
  double r = 2, v1, v2, l, result;
  while (r >= 1.0 || r < 1e-300) {
    v1 = 2.0 * randomDouble() - 1.0;
    v2 = 2.0 * randomDouble() - 1.0;
    r = v1 * v1 + v2 * v2;
  }
  l = v1 * sqrt(-2.0 * ::log(r) / r);
  result = avg + std * l;
  return (result);
}

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 48 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 11 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 "s
          << 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()

string eonc::helpers::getRelevantFile

(

std::string

filename

)

Definition at line 76 of file HelperFunctions.cpp.

                                                   {
  string filenameRelevant;
  string filenamePrefix;
  string filenamePostfix;
 
  // check if the _cp version of the file is present
  int i = filename.rfind(".");
  filenamePrefix.assign(filename, 0, i);
  filenamePostfix.assign(filename, i, filename.size());
  filenameRelevant = filenamePrefix + "_cp" + filenamePostfix;
  if (existsFile(filenameRelevant)) {
    return filenameRelevant;
  }
  // check if the _in version of the file is present
  filenameRelevant = filenamePrefix + "_in" + filenamePostfix;
  if (existsFile(filenameRelevant)) {
    return filenameRelevant;
  }
  // otherwise return original filename
  return filename;
}

getTime()

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

Definition at line 45 of file HelperFunctions.cpp.

                                                                 {
  // Wall-clock time via C++11 chrono (portable)
  using namespace std::chrono;
  auto now = steady_clock::now();
  *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
}

identical()

bool eonc::geometry::identical	(	const Matter &	m1,
		const Matter &	m2,
		const double	distanceDifference )

Definition at line 37 of file GeometryAnalysis.cpp.

                                                                {
 
  AtomMatrix r1 = m1.getPositions();
  AtomMatrix r2 = m2.getPositions();
 
  std::set<int> matched;
  double tolerance = distanceDifference;
 
  if (r1.rows() != r2.rows()) {
    return false;
  }
  int N = r1.rows();
 
  for (int i = 0; i < N; i++) {
    if (std::fabs((m1.pbc(r1.row(i) - r2.row(i))).norm()) < tolerance &&
        m1.getAtomicNr(i) == m2.getAtomicNr(i)) {
      matched.insert(i);
    }
  }
 
  for (int j = 0; j < N; j++) {
 
    if (matched.count(j) == 1)
      continue;
 
    for (int k = 0; k < N; k++) {
      if (matched.count(j) == 1)
        break;
 
      if (std::fabs((m1.pbc(r1.row(j) - r2.row(k))).norm()) < tolerance &&
          m1.getAtomicNr(j) == m2.getAtomicNr(k)) {
        matched.insert(j);
      }
    }
  }
 
  if (matched.size() == static_cast<unsigned>(N)) {
    return true;
  } else {
    return false;
  }
}

isNumber()

bool eonc::helpers::isNumber

(

std::string_view

token

)

Figure out if a string has a number in it.

Based on

Definition at line 55 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 98 of file HelperFunctions.cpp.

                                                            {
  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 121 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 136 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 151 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.
    const AtomMatrix &initPos = initial.getPositions();
    AtomMatrix pos = target.getPositionsCopy();
    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);
    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()

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

(

std::unique_ptr< Parameters >

params

)

Definition at line 39 of file Job.cpp.

                                                             {
  switch (params->main_options.job) {
    using enum JobType;
  case Process_Search: {
    return (std::make_unique<ProcessSearchJob>(std::move(params)));
    break;
  }
  case Saddle_Search: {
    return (std::make_unique<SaddleSearchJob>(std::move(params)));
    break;
  }
  case Minimization: {
    return (std::make_unique<MinimizationJob>(std::move(params)));
    break;
  }
  case Point: {
    return (std::make_unique<PointJob>(std::move(params)));
    break;
  }
  case Parallel_Replica: {
    return (std::make_unique<ParallelReplicaJob>(std::move(params)));
    break;
  }
  case Safe_Hyperdynamics: {
    return (std::make_unique<SafeHyperJob>(std::move(params)));
    break;
  }
  case TAD: {
    return (std::make_unique<TADJob>(std::move(params)));
    break;
  }
  case Replica_Exchange: {
    return (std::make_unique<ReplicaExchangeJob>(std::move(params)));
    break;
  }
  case Basin_Hopping: {
    return (std::make_unique<BasinHoppingJob>(std::move(params)));
    break;
  }
  case Hessian: {
    return (std::make_unique<HessianJob>(std::move(params)));
    break;
  }
  case Finite_Difference: {
    return (std::make_unique<FiniteDifferenceJob>(std::move(params)));
    break;
  }
  case Nudged_Elastic_Band: {
    return (std::make_unique<NudgedElasticBandJob>(std::move(params)));
    break;
  }
  case Dynamics: {
    return (std::make_unique<DynamicsJob>(std::move(params)));
    break;
  }
  case Prefactor: {
    return (std::make_unique<PrefactorJob>(std::move(params)));
    break;
  }
  case Global_Optimization: {
    return (std::make_unique<GlobalOptimizationJob>(std::move(params)));
    break;
  }
  case Structure_Comparison: {
    return (std::make_unique<StructureComparisonJob>(std::move(params)));
    break;
  }
  case Monte_Carlo: {
    return (std::make_unique<MonteCarloJob>(std::move(params)));
    break;
  }
#ifdef WITH_GP_SURROGATE
  case GP_Surrogate: {
    return (std::make_unique<GPSurrogateJob>(std::move(params)));
    break;
  }
#endif
  case OH_TST: {
    return (std::make_unique<OHTSTJob>(std::move(params)));
  }
  default:
    throw std::runtime_error("No known job could be constructed");
    break;
  }
}

makeOrthogonal()

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

Definition at line 40 of file HelperFunctions.cpp.

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

makePotential() [1/2]

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

(

const Parameters &

params

)

Definition at line 118 of file Potential.cpp.

                                                                 {
  // Inject config-file path before any potential constructor runs
  PluginLoader::instance().add_config_paths(
      params.potential_options.potentialsPath);
  return makePotential(params.potential_options.potential, params);
}

makePotential() [2/2]

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

Definition at line 124 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
  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
  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;
  }
}

maxAtomMotion()

double eonc::geometry::maxAtomMotion

(

const AtomMatrix

v1

)

Definition at line 29 of file GeometryAnalysis.cpp.

                                                      {
  return v1.rowwise().norm().maxCoeff();
}

maxAtomMotionApplied()

AtomMatrix eonc::geometry::maxAtomMotionApplied	(	const AtomMatrix	v1,
		double	maxMotion )

Definition at line 32 of file GeometryAnalysis.cpp.

                                                                  {
  AtomMatrix v2(v1);
 
  double max = maxAtomMotion(v1);
  if (max > maxMotion) {
    v2 *= maxMotion / max;
  }
  return v2;
}

maxAtomMotionAppliedV()

VectorXd eonc::geometry::maxAtomMotionAppliedV	(	const VectorXd	v1,
		double	maxMotion )

Definition at line 33 of file GeometryAnalysis.cpp.

                                                                 {
  VectorXd v2(v1);
 
  double max = maxAtomMotionV(v1);
  if (max > maxMotion) {
    v2 *= maxMotion / max;
  }
  return v2;
}

maxAtomMotionV()

double eonc::geometry::maxAtomMotionV

(

const VectorXd

v1

)

Definition at line 30 of file GeometryAnalysis.cpp.

                                                     {
  double max = 0.0;
  long n = v1.rows();
  if (n < 3) {
    // Vector too short for 3D atom grouping; treat as single displacement
    return v1.norm();
  }
  for (long i = 0; i + 3 <= n; i += 3) {
    double norm = v1.segment<3>(i).norm();
    if (max < norm) {
      max = norm;
    }
  }
  // Handle trailing elements (vector size not multiple of 3)
  long rem = n % 3;
  if (rem > 0) {
    double norm = v1.tail(rem).norm();
    if (max < norm) {
      max = norm;
    }
  }
  return max;
}

maxMotionApplied()

AtomMatrix eonc::geometry::maxMotionApplied	(	const AtomMatrix	v1,
		double	maxMotion )

Definition at line 34 of file GeometryAnalysis.cpp.

                                                              {
  AtomMatrix v2(v1);
 
  double max = v1.norm();
  if (max > maxMotion) {
    v2 *= maxMotion / max;
  }
  return v2;
}

maxMotionAppliedV()

VectorXd eonc::geometry::maxMotionAppliedV	(	const VectorXd	v1,
		double	maxMotion )

Definition at line 35 of file GeometryAnalysis.cpp.

                                                             {
  VectorXd v2(v1);
 
  double max = v1.norm();
  if (max > maxMotion) {
    v2 *= maxMotion / max;
  }
  return v2;
}

numAtomsMoved()

long eonc::geometry::numAtomsMoved	(	const AtomMatrix	v1,
		double	cutoff )

Definition at line 31 of file GeometryAnalysis.cpp.

                                                                   {
  long num = 0;
  for (int i = 0; i < v1.rows(); i++) {
    double norm = v1.row(i).norm();
    if (norm >= cutoff) {
      num += 1;
    }
  }
  return num;
}

projectOutRotTrans()

void eonc::geometry::projectOutRotTrans	(	Eigen::VectorXd &	step,
		const AtomMatrix &	positions )

Definition at line 23 of file GeometryAnalysis.cpp.

                                                                     {
  long nAtoms = positions.rows();
  long dof = nAtoms * 3;
 
  // 1. Compute center of mass (unweighted geometric center)
  Eigen::Vector3d com = Eigen::Vector3d::Zero();
  for (long i = 0; i < nAtoms; ++i) {
    com(0) += positions(i, 0);
    com(1) += positions(i, 1);
    com(2) += positions(i, 2);
  }
  com /= static_cast<double>(nAtoms);
 
  // 2. Construct 6 basis vectors for rigid-body translation and rotation
  std::vector<Eigen::VectorXd> basis;
  basis.reserve(6);
 
  // Translational basis vectors
  for (int d = 0; d < 3; ++d) {
    Eigen::VectorXd t = Eigen::VectorXd::Zero(dof);
    for (long j = 0; j < nAtoms; ++j) {
      t(3 * j + d) = 1.0;
    }
    basis.push_back(t);
  }
 
  // Rotational basis vectors (infinitesimal rotations around COM)
  Eigen::VectorXd rx = Eigen::VectorXd::Zero(dof);
  Eigen::VectorXd ry = Eigen::VectorXd::Zero(dof);
  Eigen::VectorXd rz = Eigen::VectorXd::Zero(dof);
 
  for (long i = 0; i < nAtoms; ++i) {
    double x = positions(i, 0) - com(0);
    double y = positions(i, 1) - com(1);
    double z = positions(i, 2) - com(2);
    // Rotation around x-axis: cross(xhat, r) = (0, -z, y)
    rx(3 * i + 1) = -z;
    rx(3 * i + 2) = y;
    // Rotation around y-axis: cross(yhat, r) = (z, 0, -x)
    ry(3 * i + 0) = z;
    ry(3 * i + 2) = -x;
    // Rotation around z-axis: cross(zhat, r) = (-y, x, 0)
    rz(3 * i + 0) = -y;
    rz(3 * i + 1) = x;
  }
  basis.push_back(rx);
  basis.push_back(ry);
  basis.push_back(rz);
 
  // 3. Modified Gram-Schmidt orthonormalization
  std::vector<Eigen::VectorXd> ortho;
  ortho.reserve(6);
 
  for (auto &v : basis) {
    Eigen::VectorXd u = v;
    for (const auto &e : ortho) {
      u -= u.dot(e) * e;
    }
    // Handles linear molecules where one rotational mode is degenerate
    if (u.norm() > 1e-9) {
      u.normalize();
      ortho.push_back(u);
    }
  }
 
  // 4. Project out rigid-body components from step
  for (const auto &e : ortho) {
    step -= step.dot(e) * e;
  }
}

pushApart()

void eonc::geometry::pushApart	(	std::shared_ptr< Matter >	m1,
		double	minDistance )

Definition at line 41 of file GeometryAnalysis.cpp.

                                                                         {
  if (minDistance <= 0)
    return;
 
  AtomMatrix r1 = m1->getPositions();
  AtomMatrix Force(r1.rows(), 3);
  double f = 0.025;
  double cut = minDistance;
  double pushAparts = 500;
  for (int p = 0; p < r1.rows(); p++) {
    for (int axis = 0; axis <= 2; axis++) {
      Force(p, axis) = 0;
    }
  }
  for (int count = 0; count < pushAparts; count++) {
    int moved = 0;
    for (int i = 0; i < r1.rows(); i++) {
      for (int j = i + 1; j < r1.rows(); j++) {
        double d = m1->distance(i, j);
        if (d < cut) {
          moved++;
          for (int axis = 0; axis <= 2; axis++) {
            double componant = f * (r1(i, axis) - r1(j, axis)) / d;
            Force(i, axis) += componant;
            Force(j, axis) -= componant;
          }
        }
      }
    }
    if (moved == 0)
      break;
    for (int k = 0; k < r1.rows(); k++) {
      for (int axis = 0; axis <= 2; axis++) {
        r1(k, axis) += Force(k, axis);
        Force(k, axis) = 0;
      }
    }
    m1->setPositions(r1);
  }
}

random()

double eonc::rng::random

(

long

newSeed = 0

)

Definition at line 35 of file RandomNumbers.cpp.

                                   {
  static long seed = -1;
  if (newSeed) {
    seed = -newSeed;
  }
  int j;
  long k;
  static long seed2 = 123456789;
  static long iy = 0;
  static long iv[NTAB];
  double temp;
  if (seed <= 0) {
    if (-(seed) < 1)
      seed = 3;
    else
      seed = -(seed);
    seed2 = (seed);
    for (j = NTAB + 7; j >= 0; j--) {
      k = (seed) / IQ1;
      seed = IA1 * (seed - k * IQ1) - k * IR1;
      if (seed < 0)
        seed += IM1;
      if (j < NTAB)
        iv[j] = seed;
    }
    iy = iv[0];
  }
  k = (seed) / IQ1;
  seed = IA1 * (seed - k * IQ1) - k * IR1;
  if (seed < 0)
    seed += IM1;
  k = seed2 / IQ2;
  seed2 = IA2 * (seed2 - k * IQ2) - k * IR2;
  if (seed2 < 0)
    seed2 += IM2;
  j = int(iy / NDIV);
  iy = iv[j] - seed2;
  iv[j] = seed;
  if (iy < 1)
    iy += IMM1;
  if ((temp = double(AM * iy)) > RNMX)
    return RNMX;
  else
    return temp;
}

randomDouble()

double eonc::rng::randomDouble

(

)

Definition at line 36 of file RandomNumbers.cpp.

{ return (random()); }

randomInt()

long eonc::rng::randomInt	(	int	lower,
		int	upper )

Definition at line 40 of file RandomNumbers.cpp.

                                            {
  return lround((upper - lower) * randomDouble() + lower);
}

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 314 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 263 of file HelperFunctions.cpp.

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

rotationExtract()

RotationMatrix eonc::geometry::rotationExtract	(	const AtomMatrix	r1,
		const AtomMatrix	r2 )

Definition at line 21 of file GeometryAnalysis.cpp.

                                                                    {
  RotationMatrix R;
 
  // Determine optimal rotation
  // Horn, J. Opt. Soc. Am. A, 1987
  Matrix3d m = r1.transpose() * r2;
 
  double sxx = m(0, 0);
  double sxy = m(0, 1);
  double sxz = m(0, 2);
  double syx = m(1, 0);
  double syy = m(1, 1);
  double syz = m(1, 2);
  double szx = m(2, 0);
  double szy = m(2, 1);
  double szz = m(2, 2);
 
  Matrix4d n;
  n.setZero();
  n(0, 1) = syz - szy;
  n(0, 2) = szx - sxz;
  n(0, 3) = sxy - syx;
 
  n(1, 2) = sxy + syx;
  n(1, 3) = szx + sxz;
 
  n(2, 3) = syz + szy;
 
  n += n.transpose().eval();
 
  n(0, 0) = sxx + syy + szz;
  n(1, 1) = sxx - syy - szz;
  n(2, 2) = -sxx + syy - szz;
  n(3, 3) = -sxx - syy + szz;
 
  Eigen::SelfAdjointEigenSolver<Matrix4d> es(n);
  Eigen::Vector4d maxv = es.eigenvectors().col(3);
 
  double aa = maxv[0] * maxv[0];
  double bb = maxv[1] * maxv[1];
  double cc = maxv[2] * maxv[2];
  double dd = maxv[3] * maxv[3];
  double ab = maxv[0] * maxv[1];
  double ac = maxv[0] * maxv[2];
  double ad = maxv[0] * maxv[3];
  double bc = maxv[1] * maxv[2];
  double bd = maxv[1] * maxv[3];
  double cd = maxv[2] * maxv[3];
 
  R(0, 0) = aa + bb - cc - dd;
  R(0, 1) = 2 * (bc - ad);
  R(0, 2) = 2 * (bd + ac);
  R(1, 0) = 2 * (bc + ad);
  R(1, 1) = aa - bb + cc - dd;
  R(1, 2) = 2 * (cd - ab);
  R(2, 0) = 2 * (bd - ac);
  R(2, 1) = 2 * (cd + ab);
  R(2, 2) = aa - bb - cc + dd;
 
  return R;
}

rotationMatch()

bool eonc::geometry::rotationMatch	(	const Matter &	m1,
		const Matter &	m2,
		const double	max_diff )

Definition at line 22 of file GeometryAnalysis.cpp.

                                                          {
  AtomMatrix r1 = m1.getPositions();
  AtomMatrix r2 = m2.getPositions();
 
  // Align centroids
  Eigen::VectorXd c1(3);
  Eigen::VectorXd c2(3);
 
  c1[0] = r1.col(0).sum();
  c1[1] = r1.col(1).sum();
  c1[2] = r1.col(2).sum();
  c2[0] = r2.col(0).sum();
  c2[1] = r2.col(1).sum();
  c2[2] = r2.col(2).sum();
  c1 /= r1.rows();
  c2 /= r2.rows();
 
  for (int i = 0; i < r1.rows(); i++) {
    r1(i, 0) -= c1[0];
    r1(i, 1) -= c1[1];
    r1(i, 2) -= c1[2];
 
    r2(i, 0) -= c2[0];
    r2(i, 1) -= c2[1];
    r2(i, 2) -= c2[2];
  }
 
  RotationMatrix R = rotationExtract(r1, r2);
 
  // Eigen is transposed relative to numpy
  r2 = r2 * R;
 
  for (int i = 0; i < r1.rows(); i++) {
    double diff = (r2.row(i) - r1.row(i)).norm();
    if (diff > max_diff) {
      return false;
    }
  }
  return true;
}

rotationRemove()

void eonc::geometry::rotationRemove	(	const AtomMatrix	r1,
		std::shared_ptr< Matter >	m2 )

Definition at line 24 of file GeometryAnalysis.cpp.

                                                              {
  // Skip for extended systems (slabs/surfaces with frozen atoms).
  // Rigid-body rotation and translation are not well-defined when the
  // system is anchored by frozen atoms.
  if (m2->numberOfFixedAtoms() > 0) {
    return;
  }
 
  AtomMatrix r2 = m2->getPositions();
  long n = r1_passed.size();
 
  // Compute displacement as flat 3N vector
  Eigen::VectorXd step(n);
  Eigen::Map<const Eigen::VectorXd> r1_flat(r1_passed.data(), n);
  Eigen::Map<const Eigen::VectorXd> r2_flat(r2.data(), n);
  step = r2_flat - r1_flat;
 
  // Project out rigid-body translation and rotation
  projectOutRotTrans(step, r1_passed);
 
  // Reconstruct positions: r1 + projected step
  Eigen::VectorXd result = r1_flat + step;
  AtomMatrix resultMat(r1_passed.rows(), 3);
  Eigen::Map<Eigen::VectorXd>(resultMat.data(), n) = result;
 
  m2->setPositions(resultMat);
}

saveMode() [1/2]

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

Definition at line 211 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 200 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;
}

sortedR()

bool eonc::geometry::sortedR	(	const Matter &	m1,
		const Matter &	m2,
		const double	distanceDifference )

Definition at line 39 of file GeometryAnalysis.cpp.

                                                              {
  EONC_LOG_INFO("In sortedR");
  AtomMatrix r1 = m1.getPositions();
  AtomMatrix r2 = m2.getPositions();
  double tolerance = distanceDifference;
  int matches = 0;
 
  if (r1.rows() != r2.rows()) {
    return false;
  }
 
  // Allocate memory for rdf1 and rdf2
  std::vector<std::set<atom, by_atom>> rdf1(r1.rows());
  std::vector<std::set<atom, by_atom>> rdf2(r2.rows());
 
  for (int i2 = 0; i2 < r2.rows(); i2++) {
    rdf2[i2].clear();
    for (int j2 = 0; j2 < r2.rows(); j2++) {
      if (j2 == i2)
        continue;
      atom a2;
      a2.r = m2.distance(i2, j2);
      a2.z = m2.getAtomicNr(j2);
      rdf2[i2].insert(a2);
      rdf2[j2].insert(a2);
    }
  }
 
  for (int i1 = 0; i1 < r1.rows(); i1++) {
    if (matches == i1 - 2) {
      return false;
    }
    for (int j1 = 0; j1 < r1.rows(); j1++) {
      if (j1 == i1)
        continue;
      atom a;
      a.r = m1.distance(i1, j1);
      a.z = m1.getAtomicNr(j1);
      rdf1[i1].insert(a);
      rdf1[j1].insert(a);
    }
    for (int x = 0; x < r2.rows(); x++) {
      auto it2 = rdf2[x].begin();
      auto it = rdf1[i1].begin();
      int c = 0;
      int counter = 0;
      for (; c < r1.rows(); c++) {
        if (it == rdf1[i1].end() || it2 == rdf2[x].end())
          break;
        atom k1 = *it;
        atom k2 = *it2;
        if (std::fabs(k1.r - k2.r) < tolerance && k1.z == k2.z) {
          counter++;
        } else {
          EONC_LOG_INFO("No match");
          break;
        }
        ++it;
        ++it2;
      }
      if (counter == r1.rows()) {
        matches++;
      } else {
        EONC_LOG_INFO("No match");
      }
    }
  }
 
  return matches >= r1.rows();
}

split_string_int()

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

Definition at line 224 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;
}

translationRemove()

void eonc::geometry::translationRemove	(	Matter &	m1,
		const AtomMatrix	r1 )

Definition at line 27 of file GeometryAnalysis.cpp.

                                                                           {
  AtomMatrix r1 = m1.getPositions();
  AtomMatrix r2 = r2_passed;
 
  // net displacement
  Eigen::VectorXd disp(3);
  AtomMatrix r12 = m1.pbc(r2 - r1);
 
  disp[0] = r12.col(0).sum();
  disp[1] = r12.col(1).sum();
  disp[2] = r12.col(2).sum();
  disp /= r1.rows();
 
  for (int i = 0; i < r1.rows(); i++) {
    r1(i, 0) += disp[0];
    r1(i, 1) += disp[1];
    r1(i, 2) += disp[2];
  }
 
  m1.setPositions(r1);
  return;
}

Variable Documentation

pi

double eonc::helpers::pi = 3.14159265358979323846

inlineconstexpr

Definition at line 28 of file HelperFunctions.h.