eonc::BasinHoppingJob Class Reference

#include <BasinHoppingJob.h>

Inheritance diagram for eonc::BasinHoppingJob:

Public Member Functions
	BasinHoppingJob (std::unique_ptr< Parameters > parameters, Runtime &rt)
	~BasinHoppingJob (void)=default
std::vector< std::string >	run (void) override
	Virtual run; used solely for dynamic dispatch.
Public Member Functions inherited from eonc::Job
void	adoptRuntime (std::unique_ptr< Runtime > rt)
	Take ownership of a Runtime previously passed as Runtime&.
	Job (std::unique_ptr< Parameters > parameters, Runtime &rt)
	Borrow: caller keeps Runtime alive (CLI stack / Python Session).
	Job (std::unique_ptr< Parameters > parameters, std::unique_ptr< Runtime > rt)
	Own a Runtime (one-shot makeJob / rvalue).
	Job (std::unique_ptr< Parameters > parameters)
	Own a default-constructed Runtime.
	Job (std::shared_ptr< Potential > potPassed, const Parameters &parameters)
virtual	~Job ()=default
JobType	getType ()
PotRegistry &	pots () noexcept
void	releasePotential ()
	Drop the Potential so on_destroyed is recorded before Runtime dies.

Static Public Member Functions
static double	metropolisProbability (double de, double kB, double temperature)
	Metropolis weight for energy change de.

Private Member Functions
VectorXd	calculateDistanceFromCenter (Matter *matter)
AtomMatrix	displaceRandom (double maxDisplacement)
void	randomSwap (Matter *matter)
std::vector< long >	getElements (Matter *matter)

Private Attributes
std::shared_ptr< Matter >	current
std::shared_ptr< Matter >	trial
std::vector< std::string >	returnFiles
int	jump_count {0}
int	disp_count {0}
int	swap_count {0}
int	fcalls {0}
std::vector< std::shared_ptr< Matter > >	uniqueStructures
std::vector< double >	uniqueEnergies
eonc::log::Scoped	log

Friends
class	BasinHoppingDisplaceAccess
struct	BasinHoppingElementsTest

Additional Inherited Members
Protected Attributes inherited from eonc::Job
JobType	jtype
Parameters	params
std::unique_ptr< Runtime >	owned_runtime_
	Non-null when this Job owns the composition root (one-shot makeJob).
Runtime *	runtime_
	Always valid: either owned_runtime_.get() or a caller-owned Runtime.
std::shared_ptr< Potential >	pot

Detailed Description

Definition at line 22 of file BasinHoppingJob.h.

Constructor & Destructor Documentation

BasinHoppingJob()

eonc::BasinHoppingJob::BasinHoppingJob ( std::unique_ptr< Parameters > parameters, Runtime & rt )

inline

Definition at line 24 of file BasinHoppingJob.h.

      : Job(std::move(parameters), rt),
        current{std::make_shared<Matter>(pot, params)},
        trial{std::make_shared<Matter>(pot, params)}, fcalls{0} {}

~BasinHoppingJob()

eonc::BasinHoppingJob::~BasinHoppingJob ( void )

default

Member Function Documentation

calculateDistanceFromCenter()

VectorXd eonc::BasinHoppingJob::calculateDistanceFromCenter ( Matter * matter )

private

Definition at line 475 of file BasinHoppingJob.cpp.

                                                                    {
  AtomMatrix pos = matter->getPositions();
  Vector3d cen(0, 0, 0);
  int num = matter->numberOfAtoms();
 
  cen = pos.colwise().sum() / static_cast<double>(num);
 
  VectorXd dist(num);
 
  for (int n = 0; n < num; n++) {
    pos.row(n) -= cen;
    dist(n) = pos.row(n).norm();
  }
 
  return dist;
}

displaceRandom()

AtomMatrix eonc::BasinHoppingJob::displaceRandom ( double maxDisplacement )

private

Definition at line 333 of file BasinHoppingJob.cpp.

                                                                 {
  disp_count++;
  // Create a random displacement
  AtomMatrix displacement;
  displacement.resize(trial->numberOfAtoms(), 3);
  displacement.setZero();
  VectorXd distvec = calculateDistanceFromCenter(current.get());
  // Coincident atoms share a zero radius. Dividing by it writes NaN and
  // the hop cannot leave the center. Those atoms are the outer shell, so
  // the step stays the unscaled displacement.
  const double radius = distvec.size() > 0 ? distvec.maxCoeff() : 0.0;
  int num = trial->numberOfAtoms();
  int m = 0;
  if (params.basin_hopping_options().single_atom_displace) {
    m = eonc::rng::randomInt(0, trial->numberOfAtoms() - 1);
    num = m + 1;
  }
 
  for (int i = m; i < num; i++) {
    double dist = distvec(i);
    double disp = 0.0; // displacement size, possibly scaled
 
    if (!trial->getFixed(i)) {
      const std::string &algorithm =
          params.basin_hopping_options().displacement_algorithm;
      if (algorithm == "standard") {
        disp = curDisplacement;
      }
      // scale displacement linearly with the particle radius
      else if (algorithm == "linear") {
        disp =
            radius == 0.0 ? curDisplacement : curDisplacement * (dist / radius);
      }
      // scale displacement quadratically with the particle radius
      else if (algorithm == "quadratic") {
        if (radius == 0.0) {
          disp = curDisplacement;
        } else {
          const double scale = dist / radius;
          disp = curDisplacement * scale * scale;
        }
      } else {
        log = eonc::log::traceback();
        QUILL_LOG_CRITICAL(log, "Unknown displacement_algorithm\n");
        throw std::invalid_argument(
            std::format("[Basin Hopping] unknown displacement_algorithm: {}",
                        params.basin_hopping_options().displacement_algorithm));
      }
      for (int j = 0; j < 3; j++) {
        if (params.basin_hopping_options().displacement_distribution ==
            "uniform") {
          displacement(i, j) = eonc::rng::randomDouble(2 * disp) - disp;
        } else if (params.basin_hopping_options().displacement_distribution ==
                   "gaussian") {
          displacement(i, j) = eonc::rng::gaussRandom(0.0, disp);
        } else {
          log = eonc::log::traceback();
          QUILL_LOG_CRITICAL(log, "Unknown displacement_distribution\n");
          throw std::invalid_argument(std::format(
              "[Basin Hopping] unknown displacement_distribution: {}",
              params.basin_hopping_options().displacement_distribution));
        }
      }
    }
  }
  return displacement;
}

getElements()

std::vector< long > eonc::BasinHoppingJob::getElements ( Matter * matter )

private

Definition at line 451 of file BasinHoppingJob.cpp.

                                                           {
  // Z is 0..118. A 118-slot table stores 118 one past the end.
  constexpr int kElementSlots = 119;
  std::array<int, kElementSlots> allElements{};
  std::vector<long> elements;
 
  for (long y = 0; y < matter->numberOfAtoms(); ++y) {
    if (!matter->getFixed(y)) {
      const long z = matter->getAtomicNr(y);
      if (z >= 0 && z < kElementSlots) {
        allElements[static_cast<size_t>(z)] = 1;
      }
    }
  }
 
  for (int i = 0; i < kElementSlots; ++i) {
    if (allElements[static_cast<size_t>(i)] != 0) {
      elements.push_back(i);
    }
  }
 
  return elements;
}

metropolisProbability()

double eonc::BasinHoppingJob::metropolisProbability ( double de, double kB, double temperature )

static

Metropolis weight for energy change de.

Divides by kB * temperature only for an uphill hop when both are positive. Non-positive temperature or kB refuses that hop.

Definition at line 322 of file BasinHoppingJob.cpp.

                                                                  {
  if (!(de > 0.0)) {
    return 1.0;
  }
  if (!(temperature > 0.0) || !(kB > 0.0)) {
    return 0.0;
  }
  return std::exp(-de / (kB * temperature));
}

randomSwap()

void eonc::BasinHoppingJob::randomSwap ( Matter * matter )

private

Definition at line 401 of file BasinHoppingJob.cpp.

                                               {
  swap_count++;
  std::vector<long> Elements;
  Elements = getElements(matter);
 
  long ela;
  long elb;
  long ia = eonc::rng::randomInt(0, Elements.size() - 1);
  ela = Elements.at(ia);
  Elements.erase(Elements.begin() + ia);
 
  long ib = eonc::rng::randomInt(0, Elements.size() - 1);
  elb = Elements.at(ib);
 
  int changera = 0;
  int changerb = 0;
 
  changera = eonc::rng::randomInt(0, matter->numberOfAtoms() - 1);
  int guard = 0;
  while ((matter->getAtomicNr(changera) != ela || matter->getFixed(changera)) &&
         guard < 10000) {
    changera = eonc::rng::randomInt(0, matter->numberOfAtoms() - 1);
    guard++;
  }
  changerb = eonc::rng::randomInt(0, matter->numberOfAtoms() - 1);
  guard = 0;
  while ((matter->getAtomicNr(changerb) != elb || matter->getFixed(changerb) ||
          changerb == changera) &&
         guard < 10000) {
    changerb = eonc::rng::randomInt(0, matter->numberOfAtoms() - 1);
    guard++;
  }
  if (matter->getAtomicNr(changera) != ela || matter->getFixed(changera) ||
      matter->getAtomicNr(changerb) != elb || matter->getFixed(changerb)) {
    return;
  }
 
  double posax = matter->getPosition(changera, 0);
  double posay = matter->getPosition(changera, 1);
  double posaz = matter->getPosition(changera, 2);
 
  matter->setPosition(changera, 0, matter->getPosition(changerb, 0));
  matter->setPosition(changera, 1, matter->getPosition(changerb, 1));
  matter->setPosition(changera, 2, matter->getPosition(changerb, 2));
 
  matter->setPosition(changerb, 0, posax);
  matter->setPosition(changerb, 1, posay);
  matter->setPosition(changerb, 2, posaz);
}

run()

std::vector< std::string > eonc::BasinHoppingJob::run ( void )

overridevirtual

Virtual run; used solely for dynamic dispatch.

Implements eonc::Job.

Definition at line 32 of file BasinHoppingJob.cpp.

                                          {
  bool swapMove;
  double swap_accept = 0.0;
  jump_count = 0; // count of jump movies
  swap_count = 0; // count of swap moves
  disp_count = 0; // count of displacement moves
  // Quench tail may not run when stop_energy breaks first.
  int quench_displacements = 0;
  int consecutive_rejected_trials = 0;
  double totalAccept = 0.0;
  std::unique_ptr<Matter> minTrial = std::make_unique<Matter>(pot, params);
  std::unique_ptr<Matter> swapTrial = std::make_unique<Matter>(pot, params);
 
  std::string conFilename =
      eonc::helpers::getRelevantFile(params.main_options().conFilename);
  if (!eonc::io::io_ok(current->con2matter(conFilename))) {
    QUILL_LOG_CRITICAL(log, "Failed to load {}", conFilename);
    throw std::runtime_error("failed to load " + conFilename);
  }
 
  // Sanity Check
  std::vector<long> Elements;
  Elements = getElements(current.get());
  if (params.basin_hopping_options().swap_probability > 0 &&
      Elements.size() == 1) {
    log = eonc::log::traceback();
    QUILL_LOG_CRITICAL(log,
                       "error: [Basin Hopping] swap move probability must be "
                       "zero if there is only one element type\n");
    throw std::invalid_argument(
        "[Basin Hopping] swap_probability must be zero with one element type");
  }
 
  double randomProb =
      params.basin_hopping_options().initial_random_structure_probability;
  if (randomProb > 0.0) {
    QUILL_LOG_DEBUG(log, "generating random structure with probability {:.4f}",
                    randomProb);
  }
  double u = eonc::rng::random();
  if (u < params.basin_hopping_options().initial_random_structure_probability) {
    AtomMatrix randomPositions = current->getPositionsFree();
    for (int i = 0; i < current->numberOfFreeAtoms(); i++) {
      for (int j = 0; j < 3; j++) {
        randomPositions(i, j) = eonc::rng::random();
      }
    }
    randomPositions *= current->getCell();
    current->setPositionsFree(randomPositions);
 
    eonc::geometry::pushApart(
        current, params.basin_hopping_options().push_apart_distance);
  }
 
  *trial = *current;
  *minTrial = *current;
 
  current->relax(true);
 
  double currentEnergy = current->getPotentialEnergy();
  double minimumEnergy = currentEnergy;
 
  auto minimumEnergyStructure = std::make_shared<Matter>(pot, params);
  *minimumEnergyStructure = *current;
  int nsteps = params.basin_hopping_options().steps +
               params.basin_hopping_options().quenching_steps;
 
  QUILL_LOG_DEBUG(
      log, "[Basin Hopping] {:4s} {:12s} {:12s} {:12s} {:4s} {:5s} {:5s}",
      "step", "current", "trial", "global min", "fc", "ar", "md");
  QUILL_LOG_DEBUG(
      log, "[Basin Hopping] {:4s} {:12s} {:12s} {:12s} {:4s} {:5s} {:5s}",
      "----", "-------", "-----", "----------", "--", "--", "--");
 
  int recentAccept = 0;
  double curDisplacement = params.basin_hopping_options().displacement;
 
  for (int step = 0; step < nsteps; step++) {
 
    // Swap or displace
    if (eonc::rng::randomDouble(1.0) <
            params.basin_hopping_options().swap_probability &&
        step < params.basin_hopping_options().steps) {
      *swapTrial = *current;
      randomSwap(swapTrial.get());
      swapMove = true;
      *minTrial = *swapTrial;
    } else {
      AtomMatrix displacement;
      displacement = displaceRandom(curDisplacement);
      if (step >= params.basin_hopping_options().steps) {
        quench_displacements++;
      }
 
      trial->setPositions(current->getPositions() + displacement);
      swapMove = false;
      eonc::geometry::pushApart(
          trial, params.basin_hopping_options().push_apart_distance);
 
      *minTrial = *trial;
    }
 
    if (params.debug_options().write_movies) {
      if (!eonc::io::io_ok(trial->matter2con("trials", true))) {
        QUILL_LOG_WARNING(log, "Failed to append trials movie frame");
      }
    }
 
    minTrial->relax(true);
 
    double deltaE = minTrial->getPotentialEnergy() - currentEnergy;
    double p = 0.0;
    if (step >= params.basin_hopping_options().steps) {
      if (deltaE <= 0.0) {
        p = 1.0;
      }
    } else {
      // Divide only for an uphill hop at positive temperature.
      p = metropolisProbability(deltaE, params.constants().kB,
                                params.main_options().temperature);
    }
 
    bool accepted = false;
    if (eonc::rng::randomDouble(1.0) < p) {
      accepted = true;
      if (params.basin_hopping_options().significant_structure) {
        *current = *minTrial;
      } else if (swapMove) {
        *current = *swapTrial;
      } else {
        *current = *trial;
      }
      if (swapMove) {
        swap_accept += 1;
      }
      if (step < params.basin_hopping_options().steps) {
        totalAccept += 1;
        recentAccept += 1;
      }
 
      currentEnergy = minTrial->getPotentialEnergy();
 
      if (currentEnergy < minimumEnergy) {
        minimumEnergy = currentEnergy;
        *minimumEnergyStructure = *minTrial;
        if (!eonc::io::io_ok(minimumEnergyStructure->matter2con("min.con"))) {
          QUILL_LOG_WARNING(log, "Failed to write min.con");
        }
      }
 
      if (params.basin_hopping_options().write_unique) {
        bool newStructure = true;
        for (unsigned int i = 0; i < uniqueEnergies.size(); i++) {
          // if minTrial has a different energy or a different structure
          // it is new, otherwise it is old
          if (std::fabs(currentEnergy - uniqueEnergies[i]) <
              params.structure_comparison_options().energy_difference) {
            Matter probe = *current;
            if (probe.compare(*uniqueStructures[i],
                              params.structure_comparison_options()
                                  .indistinguishable_atoms)) {
              newStructure = false;
            }
          }
        }
 
        if (newStructure) {
          uniqueEnergies.push_back(currentEnergy);
          auto currentCopy = std::make_shared<Matter>(pot, params);
          *currentCopy = *current;
          uniqueStructures.push_back(currentCopy);
 
          char fname[128];
          snprintf(fname, 128, "min_%.5i.con", step + 1);
          if (!eonc::io::io_ok(current->matter2con(fname))) {
            QUILL_LOG_WARNING(log, "Failed to write {}", fname);
          }
          returnFiles.push_back(fname);
 
          snprintf(fname, 128, "energy_%.5i.dat", step + 1);
          returnFiles.push_back(fname);
          {
            std::ofstream fh(fname);
            if (fh)
              fh << std::format("{:.10e}\n", currentEnergy);
          }
        }
      }
 
      consecutive_rejected_trials = 0; // STC: I think this should go here.
    } else {
      consecutive_rejected_trials++;
    }
 
    if (params.debug_options().write_movies) {
      if (!eonc::io::io_ok(minTrial->matter2con("movie", true))) {
        QUILL_LOG_WARNING(log, "Failed to append basin-hopping movie frame");
      }
    }
 
    if (minimumEnergy < params.basin_hopping_options().stop_energy) {
      break;
    }
 
    if (consecutive_rejected_trials ==
            params.basin_hopping_options().jump_max &&
        step < params.basin_hopping_options().steps) {
      consecutive_rejected_trials = 0;
      AtomMatrix jump;
      for (int j = 0; j < params.basin_hopping_options().jump_steps; j++) {
        jump_count++;
        jump = displaceRandom(curDisplacement);
        current->setPositions(current->getPositions() + jump);
        // Only a minimized jump is a basin. The raw geometry must not become
        // the Metropolis reference or the stored global minimum.
        if (params.basin_hopping_options().significant_structure) {
          eonc::geometry::pushApart(
              current, params.basin_hopping_options().push_apart_distance);
          current->relax(true);
          currentEnergy = current->getPotentialEnergy();
          if (currentEnergy < minimumEnergy) {
            minimumEnergy = currentEnergy;
            *minimumEnergyStructure = *current;
          }
        }
      }
    }
 
    int nadjust = params.basin_hopping_options().adjust_period;
    double adjustFraction = params.basin_hopping_options().adjust_fraction;
    // A zero period is not an interval; the modulo would divide by zero.
    if (nadjust != 0 && (step + 1) % nadjust == 0 &&
        params.basin_hopping_options().adjust_displacement) {
      double recentRatio =
          static_cast<double>(recentAccept) / static_cast<double>(nadjust);
      if (recentRatio > params.basin_hopping_options().target_ratio) {
        curDisplacement *= 1.0 + adjustFraction;
      } else {
        curDisplacement *= 1.0 - adjustFraction;
      }
 
      recentAccept = 0;
    }
  }
  /* Save Results */
 
  std::string resultsFilename("results.dat");
 
  if (params.debug_options().write_movies) {
    std::string movieFilename("movie.con");
    returnFiles.push_back(movieFilename);
  }
 
  {
    auto env = JobResultEnvelope::fromMinimization(
        RunStatus::GOOD, params.potential_options().potential,
        PotRegistry::get().total_force_calls(), true, minimumEnergy);
    env.job_type = "basin_hopping";
    env.random_seed = params.main_options().randomSeed;
    env.extras.emplace_back("minimum_energy", minimumEnergy);
    const double nsteps_ratio = params.basin_hopping_options().steps;
    env.extras.emplace_back("acceptance_ratio",
                            nsteps_ratio ? totalAccept / nsteps_ratio : 0.0);
    if (params.basin_hopping_options().swap_probability > 0) {
      env.extras.emplace_back(
          "swap_acceptance_ratio",
          swap_count ? swap_accept / static_cast<double>(swap_count) : 0.0);
    }
    env.extras.emplace_back(
        "total_normal_displacement_steps",
        static_cast<double>(disp_count - jump_count - quench_displacements));
    env.extras.emplace_back("total_jump_steps",
                            static_cast<double>(jump_count));
    env.extras.emplace_back("total_swap_steps",
                            static_cast<double>(swap_count));
    env.writeResultsDat(resultsFilename);
    returnFiles.push_back(resultsFilename);
  }
 
  std::string productFilename("min.con");
  if (eonc::io::io_ok(minimumEnergyStructure->matter2con(productFilename))) {
    returnFiles.push_back(productFilename);
  } else {
    QUILL_LOG_ERROR(log, "Failed to write {}", productFilename);
  }
 
  // minTrial and swapTrial automatically cleaned up by unique_ptr
  return returnFiles;
}

Friends And Related Symbol Documentation

BasinHoppingDisplaceAccess

friend class BasinHoppingDisplaceAccess

friend

Definition at line 37 of file BasinHoppingJob.h.

BasinHoppingElementsTest

friend struct BasinHoppingElementsTest

friend

Definition at line 44 of file BasinHoppingJob.h.

Member Data Documentation

current

std::shared_ptr<Matter> eonc::BasinHoppingJob::current

private

Definition at line 41 of file BasinHoppingJob.h.

disp_count

int eonc::BasinHoppingJob::disp_count {0}

private

Definition at line 47 of file BasinHoppingJob.h.

{0}; // count of displacement moves

fcalls

int eonc::BasinHoppingJob::fcalls {0}

private

Definition at line 49 of file BasinHoppingJob.h.

{0};

jump_count

int eonc::BasinHoppingJob::jump_count {0}

private

Definition at line 46 of file BasinHoppingJob.h.

{0}; // count of jump moves

log

eonc::log::Scoped eonc::BasinHoppingJob::log

private

Definition at line 53 of file BasinHoppingJob.h.

returnFiles

std::vector<std::string> eonc::BasinHoppingJob::returnFiles

private

Definition at line 45 of file BasinHoppingJob.h.

swap_count

int eonc::BasinHoppingJob::swap_count {0}

private

Definition at line 48 of file BasinHoppingJob.h.

{0}; // count of swap moves

trial

std::shared_ptr<Matter> eonc::BasinHoppingJob::trial

private

Definition at line 42 of file BasinHoppingJob.h.

uniqueEnergies

std::vector<double> eonc::BasinHoppingJob::uniqueEnergies

private

Definition at line 52 of file BasinHoppingJob.h.

uniqueStructures

std::vector<std::shared_ptr<Matter> > eonc::BasinHoppingJob::uniqueStructures

private

Definition at line 51 of file BasinHoppingJob.h.

---

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

• /home/runner/work/eOn/eOn/include/eon/BasinHoppingJob.h
• /home/runner/work/eOn/eOn/client/BasinHoppingJob.cpp