BasinHoppingJob.cpp

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Copyright (c) 2010--present, eOn Development Team
** All rights reserved.
**
** Repo:
** https://github.com/TheochemUI/eOn
*/
 
#include "eon/EonLogger.h"
#include <cmath>
#include <format>
#include <fstream>
#include <stdexcept>
#include <string>
 
#include "eon/BaseStructures.h"
#include "eon/BasinHoppingJob.h"
#include "eon/Dynamics.h"
#include "eon/HelperFunctions.h"
#include "eon/JobResult.h"
#include "eon/ObjectiveFunction.h"
#include "eon/Optimizer.h"
#include "eon/PotRegistry.h"
#include "eon/Potential.h"
 
namespace eonc {
 
std::vector<std::string> BasinHoppingJob::run() {
  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;
}
 
double BasinHoppingJob::metropolisProbability(double de, double kB,
                                              double temperature) {
  if (!(de > 0.0)) {
    return 1.0;
  }
  if (!(temperature > 0.0) || !(kB > 0.0)) {
    return 0.0;
  }
  return std::exp(-de / (kB * temperature));
}
 
AtomMatrix BasinHoppingJob::displaceRandom(double curDisplacement) {
  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;
}
 
void BasinHoppingJob::randomSwap(Matter *matter) {
  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);
}
 
std::vector<long> BasinHoppingJob::getElements(Matter *matter) {
  // 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;
}
 
VectorXd BasinHoppingJob::calculateDistanceFromCenter(Matter *matter) {
  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;
}
 
} // namespace eonc