eonc::BasinHoppingJob Class Reference

#include <BasinHoppingJob.h>

Inheritance diagram for eonc::BasinHoppingJob:

Public Member Functions
	BasinHoppingJob (std::unique_ptr< Parameters > parameters)
	~BasinHoppingJob (void)=default
std::vector< std::string >	run (void) override
	Virtual run; used solely for dynamic dispatch.
Public Member Functions inherited from eonc::Job
	Job (std::unique_ptr< Parameters > parameters)
	Job (std::shared_ptr< Potential > potPassed, const Parameters &parameters)
virtual	~Job ()=default
JobType	getType ()

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

Additional Inherited Members
Protected Attributes inherited from eonc::Job
JobType	jtype
Parameters	params
std::shared_ptr< Potential >	pot

Detailed Description

Definition at line 20 of file BasinHoppingJob.h.

Constructor & Destructor Documentation

BasinHoppingJob()

eonc::BasinHoppingJob::BasinHoppingJob ( std::unique_ptr< Parameters > parameters )

inline

Definition at line 22 of file BasinHoppingJob.h.

      : Job(std::move(parameters)),
        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 BasinHoppingJob::calculateDistanceFromCenter ( Matter * matter )

private

Definition at line 454 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 BasinHoppingJob::displaceRandom ( double maxDisplacement )

private

Definition at line 334 of file BasinHoppingJob.cpp.

                                                                 {
  disp_count++;
  // Create a random displacement
  AtomMatrix displacement;
  displacement.resize(trial->numberOfAtoms(), 3);
  displacement.setZero();
  VectorXd distvec = calculateDistanceFromCenter(current.get());
  int num = trial->numberOfAtoms();
  int m = 0;
  if (params.basin_hopping_options.single_atom_displace) {
    m = 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)) {
      if (params.basin_hopping_options.displacement_algorithm == "standard") {
        disp = curDisplacement;
      }
      // scale displacement linearly with the particle radius
      else if (params.basin_hopping_options.displacement_algorithm ==
               "linear") {
        double Cs = curDisplacement / distvec.maxCoeff();
        disp = Cs * dist;
      }
      // scale displacement quadratically with the particle radius
      else if (params.basin_hopping_options.displacement_algorithm ==
               "quadratic") {
        double Cq = curDisplacement / (distvec.maxCoeff() * distvec.maxCoeff());
        disp = Cq * dist * dist;
      } 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) = randomDouble(2 * disp) - disp;
        } else if (params.basin_hopping_options.displacement_distribution ==
                   "gaussian") {
          displacement(i, j) = 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 > BasinHoppingJob::getElements ( Matter * matter )

private

Definition at line 434 of file BasinHoppingJob.cpp.

                                                           {
  std::array<int, 118> allElements{};
  std::vector<long> elements;
 
  for (long y = 0; y < matter->numberOfAtoms(); ++y) {
    if (!matter->getFixed(y)) {
      const int index = matter->getAtomicNr(y);
      allElements[index] = 1;
    }
  }
 
  for (int i = 0; i < 118; ++i) {
    if (allElements[i] != 0) {
      elements.push_back(i);
    }
  }
 
  return elements;
}

randomSwap()

void BasinHoppingJob::randomSwap ( Matter * matter )

private

Definition at line 394 of file BasinHoppingJob.cpp.

                                               {
  swap_count++;
  std::vector<long> Elements;
  Elements = getElements(matter);
 
  long ela;
  long elb;
  long ia = randomInt(0, Elements.size() - 1);
  ela = Elements.at(ia);
  Elements.erase(Elements.begin() + ia);
 
  long ib = randomInt(0, Elements.size() - 1);
  elb = Elements.at(ib);
 
  int changera = 0;
  int changerb = 0;
 
  changera = randomInt(0, matter->numberOfAtoms() - 1);
  while (matter->getAtomicNr(changera) != ela) {
    changera = randomInt(0, matter->numberOfAtoms() - 1);
  }
 
  changerb = randomInt(0, matter->numberOfAtoms() - 1);
  while (matter->getAtomicNr(changerb) != elb) {
    changerb = randomInt(0, matter->numberOfAtoms() - 1);
  }
 
  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 > BasinHoppingJob::run ( void )

overridevirtual

Virtual run; used solely for dynamic dispatch.

Implements eonc::Job.

Definition at line 29 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
  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 = 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::helpers::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::helpers::random();
      }
    }
    randomPositions *= current->getCell();
    current->setPositionsFree(randomPositions);
 
    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;
  long totalfc;
 
  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 (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);
 
      trial->setPositions(current->getPositions() + displacement);
      swapMove = false;
      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");
      }
    }
 
    // Potential::fcalls = 0;
    minTrial->relax(true);
    // int minfcalls = Potential::fcalls;
 
    double deltaE = minTrial->getPotentialEnergy() - currentEnergy;
    double p = 0.0;
    if (step >= params.basin_hopping_options.steps) {
      if (deltaE <= 0.0) {
        p = 1.0;
      }
    } else {
      if (deltaE <= 0.0) {
        p = 1.0;
      } else {
        p = std::exp(-deltaE /
                     (params.main_options.temperature * 8.6173324e-5));
      }
    }
 
    bool accepted = false;
    if (randomDouble(1.0) < p) {
      accepted = true;
      if (params.basin_hopping_options.significant_structure) {
        *current = *minTrial;
      } 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) {
            if (current->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");
      }
    }
 
    // totalfc = Potential::fcallsTotal;
    char acceptReject[2];
    acceptReject[1] = '\0';
    if (accepted) {
      acceptReject[0] = 'A';
    } else {
      acceptReject[0] = 'R';
    }
    // QUILL_LOG_DEBUG(log, "[Basin Hopping] %5i %12.3f %12.3f %12.3f %4i
    // %5.3f %5.3f %1s\n",
    //        step+1, currentEnergy, minTrial->getPotentialEnergy(),
    //        minimumEnergy, minfcalls, totalAccept/((double)step+1),
    //        curDisplacement, acceptReject);
    // fprintf(pFile, "%6i %9ld %12.4e %12.4e\n",step+1,totalfc,currentEnergy,
    // minTrial->getPotentialEnergy());
 
    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);
        if (params.basin_hopping_options.significant_structure) {
          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;
    if ((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;
      }
 
      // QUILL_LOG_DEBUG(log, "recentRatio %.3f md: %.3f\n", recentRatio,
      // curDisplacement);
      recentAccept = 0;
    }
  }
  /* Save Results */
 
  std::string resultsFilename("results.dat");
 
  if (params.debug_options.write_movies) {
    std::string movieFilename("movie.con");
    returnFiles.push_back(movieFilename);
  }
 
  {
    std::ofstream out(resultsFilename, std::ios::binary);
    if (!out) {
      QUILL_LOG_CRITICAL(log, "Failed to open {}", resultsFilename);
      throw std::runtime_error("failed to open " + resultsFilename);
    }
    out << std::format("{} termination_reason\n", 0);
    out << "GOOD termination_reason_text\n";
    out << "basin_hopping job_type\n";
    out << std::format("{:.12e} minimum_energy\n", minimumEnergy);
    out << std::format("{} random_seed\n", params.main_options.randomSeed);
    const double nsteps_ratio = params.basin_hopping_options.steps;
    out << std::format("{:.3f} acceptance_ratio\n",
                       nsteps_ratio ? totalAccept / nsteps_ratio : 0.0);
    if (params.basin_hopping_options.swap_probability > 0) {
      out << std::format(
          "{:.3f} swap_acceptance_ratio\n",
          swap_count ? swap_accept / static_cast<double>(swap_count) : 0.0);
    }
    out << std::format("{} total_normal_displacement_steps\n",
                       disp_count - jump_count -
                           params.basin_hopping_options.quenching_steps);
    out << std::format("{} total_jump_steps\n", jump_count);
    out << std::format("{} total_swap_steps\n", swap_count);
    out << std::format("{} total_force_calls\n",
                       PotRegistry::get().total_force_calls());
    out.close();
    if (!out) {
      QUILL_LOG_CRITICAL(log, "Failed to write {}", resultsFilename);
      throw std::runtime_error("failed to write " + 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;
}

Member Data Documentation

current

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

private

Definition at line 35 of file BasinHoppingJob.h.

disp_count

int eonc::BasinHoppingJob::disp_count {0}

private

Definition at line 40 of file BasinHoppingJob.h.

{0}; // count of displacement moves

fcalls

int eonc::BasinHoppingJob::fcalls {0}

private

Definition at line 42 of file BasinHoppingJob.h.

{0};

jump_count

int eonc::BasinHoppingJob::jump_count {0}

private

Definition at line 39 of file BasinHoppingJob.h.

{0}; // count of jump moves

log

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

private

Definition at line 46 of file BasinHoppingJob.h.

returnFiles

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

private

Definition at line 38 of file BasinHoppingJob.h.

swap_count

int eonc::BasinHoppingJob::swap_count {0}

private

Definition at line 41 of file BasinHoppingJob.h.

{0}; // count of swap moves

trial

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

private

Definition at line 36 of file BasinHoppingJob.h.

uniqueEnergies

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

private

Definition at line 45 of file BasinHoppingJob.h.

uniqueStructures

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

private

Definition at line 44 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