GlobalOptimizationJob.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/GlobalOptimizationJob.h"
#include "eon/EonLogger.h"
 
#include "eon/Dynamics.h"
#include "eon/HelperFunctions.h"
#include <format>
#include <memory>
#include <stdexcept>
 
#include <algorithm>
#include <cmath>
 
namespace eonc {
 
std::vector<std::string> GlobalOptimizationJob::run() {
  std::string reactant_passed =
      eonc::helpers::getRelevantFile(params.main_options().conFilename);
  std::vector<std::string> returnFiles;
  auto matter_cur = std::make_unique<Matter>(pot, params);
  auto matter_hop = std::make_unique<Matter>(pot, params);
  if (!eonc::io::io_ok(matter_cur->con2matter(reactant_passed))) {
    QUILL_LOG_CRITICAL(log, "Failed to load {}", reactant_passed);
    throw std::runtime_error("failed to load " + reactant_passed);
  }
  bool converged;
  long nstep = params.global_optimization_options().steps;
  AtomMatrix rat_t(matter_cur->numberOfAtoms(), 3);
  QUILL_LOG_DEBUG(log, "\nBeginning minima hopping of {}",
                  reactant_passed.c_str());
  QUILL_LOG_TRACE_L1(log, "fcalls= {}", matter_cur->getForceCalls());
  QUILL_LOG_TRACE_L1(log, "epot= {:24.15E}", matter_cur->getPotentialEnergy());
  converged =
      matter_cur->relax(false, params.debug_options().write_movies,
                        params.main_options().checkpoint, "min", "matter_cur");
  QUILL_LOG_DEBUG(log, "converged {}", (converged) ? "TRUE" : "FALSE");
  earr.push_back(matter_cur->getPotentialEnergy());
  *matter_hop = *matter_cur;
  GlobalOptimizationJob::report(*matter_hop);
  for (long istep = 1; istep <= nstep; istep++) {
    // hoppingStep attempts to hop into a new state followed by a minimization
    hoppingStep(istep, *matter_cur, *matter_hop);
    // decisionStep decides to accept or reject this step if it escaped
    decisionStep(*matter_cur, *matter_hop);
    // analyzing what happened in this step if it escaped
    analyze(*matter_cur, *matter_hop);
    // reporting useful information about this hop
    report(*matter_hop);
    if (matter_cur->getPotentialEnergy() <
        params.global_optimization_options().target_energy)
      break;
  }
  for (size_t i = 0; i < earr.size(); i++) {
    double earrim1;
    if (i == 0) {
      earrim1 = earr[0];
    } else {
      earrim1 = earr[i - 1];
    }
    earrfile << std::format("{:5}  {:15.5f}  {:15.5f}  {:15.5f}  ", i + 1,
                            earr[i], earr[i] - earr[0], earr[i] - earrim1);
  }
  return returnFiles;
} // end of GlobalOptimizationJob::run
 
void GlobalOptimizationJob::analyze(Matter &matter_cur, Matter &matter_hop) {
  if (escapeResult == "failure") {
    hoppingResult = "same";
    return;
  }
  double epot = matter_hop.getPotentialEnergy();
  size_t jlo = hunt(epot);
  QUILL_LOG_TRACE_L1(log, "REZA: {}", jlo);
  if (std::abs(epot - earr[jlo]) <
      params.structure_comparison_options().energy_difference) {
    hoppingResult = "already_visited";
  } else {
    hoppingResult = "new";
  }
  if (decisionResult == "accepted") {
    double epot_hop = matter_hop.getPotentialEnergy();
    matter_cur = matter_hop;
    size_t jlo = hunt(matter_hop.getPotentialEnergy());
    if (hoppingResult == "new" && jlo == 0 && epot_hop < earr[0]) {
      QUILL_LOG_DEBUG(log,
                      "new lowest: nlmin, epot_hop, dE {:7d} {:15.5f} {:10.5f}",
                      1, epot_hop, epot_hop - earr[0]);
    }
    insert(matter_cur);
  } else if (decisionResult == "rejected") {
  } else {
    log = eonc::log::traceback();
    QUILL_LOG_CRITICAL(
        log,
        "ERROR: new minimum is neither accepted nor rejected: client stops.");
    throw std::runtime_error(
        "[Global Optimization] new minimum is neither accepted nor rejected");
  }
}
 
void GlobalOptimizationJob::examineEscape(Matter &matter_cur,
                                          Matter &matter_hop) {
  double epot, epot_hop;
  epot = matter_cur.getPotentialEnergy();
  epot_hop = matter_hop.getPotentialEnergy();
  if (std::abs(epot_hop - epot) <
      params.structure_comparison_options().energy_difference) {
    escapeResult = "failure";
  } else {
    escapeResult = "success";
  }
}
 
void GlobalOptimizationJob::applyMoveFeedbackMD() {
  if (firstStep) {
    firstStep = false;
    return;
  }
  if (hoppingResult == "same") {
    ekin *= beta1;
  } else if (hoppingResult == "already_visited") {
    ekin *= beta2;
  } else if (hoppingResult == "new") {
    ekin *= beta3;
  } else {
    log = eonc::log::traceback();
    QUILL_LOG_CRITICAL(log,
                       "ERROR: client does not know what to do with ekin.");
    QUILL_LOG_CRITICAL(log, "ERROR: client stops in applyMoveFeedbackMD.");
    throw std::runtime_error(std::format(
        "[Global Optimization] unknown hoppingResult: {}", hoppingResult));
  }
}
 
void GlobalOptimizationJob::applyDecisionFeedback() {
  if (decisionResult == "accepted") {
    ediff *= alpha1;
  } else {
    ediff *= alpha2;
  }
}
 
void GlobalOptimizationJob::report(Matter &matter_hop) {
  char C1, C2;
  double ekin_p;
  if (hoppingResult == "same") {
    C1 = 'S';
    C2 = '-';
    ekin_p = ekin * beta1;
  } else if (hoppingResult == "already_visited") {
    C1 = 'O';
    ekin_p = ekin * beta2;
  } else if (hoppingResult == "new") {
    C1 = 'N';
    ekin_p = ekin * beta3;
  } else {
    C1 = '-';
    ekin_p = ekin;
  }
  if (decisionResult == "accepted") {
    C2 = 'A';
  } else if (decisionResult == "rejected") {
    C2 = 'R';
  } else {
    C2 = '-';
  }
  double epot_hop = matter_hop.getPotentialEnergy();
  double temp = (2.0 * ekin_p / params.constants().kB);
  double dt = params.dynamics_options().time_step;
  monfile << std::format(
      "{:15.5f}  {:15.5f}  {:11}  {:12.2f}         {}{}  {:5}  {:5}", epot_hop,
      ediff, static_cast<size_t>(temp), dt, C1, C2, fcallsMove, fcallsRelax);
}
 
void GlobalOptimizationJob::decisionStep(Matter &matter_cur,
                                         Matter &matter_hop) {
  decisionResult = "unknown";
  examineEscape(matter_cur, matter_hop);
  if (escapeResult == "failure") {
    return;
  }
  if (params.global_optimization_options().decision_method == "npew") {
    acceptRejectNPEW(matter_cur, matter_hop);
  } else if (params.global_optimization_options().decision_method ==
             "boltzmann") {
    acceptRejectBoltzmann(matter_cur, matter_hop);
  } else {
    log = eonc::log::traceback();
    QUILL_LOG_CRITICAL(
        log, "ERROR: accept/reject method not specified. client stops.");
    throw std::invalid_argument(
        std::format("[Global Optimization] unknown decision_method: {}",
                    params.global_optimization_options().decision_method));
  }
  applyDecisionFeedback();
}
 
// NPEW: non-probablistic energy window (used in minima hopping method)
void GlobalOptimizationJob::acceptRejectNPEW(Matter &matter_cur,
                                             Matter &matter_hop) {
  if (matter_hop.getPotentialEnergy() <
      matter_cur.getPotentialEnergy() + ediff) {
    decisionResult = "accepted";
  } else {
    decisionResult = "rejected";
  }
}
 
void GlobalOptimizationJob::acceptRejectBoltzmann(Matter &matter_cur,
                                                  Matter &matter_hop) {
  double eTrial = matter_hop.getPotentialEnergy();
  double eCurrent = matter_cur.getPotentialEnergy();
 
  double deltaE = eTrial - eCurrent;
  const double kB = params.constants().kB;
  const double T = params.main_options().temperature;
 
  double p;
  if (deltaE <= 0.0) {
    p = 1.0;
  } else if (!(T > 0.0) || !(kB > 0.0)) {
    p = 0.0;
  } else {
    p = std::exp(-deltaE / (kB * T));
  }
 
  if (eonc::rng::randomDouble(1.0) < p) {
    decisionResult = "accepted";
  } else {
    decisionResult = "rejected";
  }
}
 
void GlobalOptimizationJob::hoppingStep(long istep, Matter &matter_cur,
                                        Matter &matter_hop) {
  bool converged;
  matter_hop = matter_cur;
  long fcalls1 = matter_hop.getForceCalls();
  if (params.global_optimization_options().move_method == "md") {
    applyMoveFeedbackMD();
    mdescape(matter_hop);
  } else if (params.global_optimization_options().move_method == "random") {
    randomMove(matter_hop);
  }
  long fcalls2 = matter_hop.getForceCalls();
  hoppingResult = "unknown";
  converged =
      matter_hop.relax(true, params.debug_options().write_movies,
                       params.main_options().checkpoint, "min", "matter_hop");
  QUILL_LOG_DEBUG(log, "converged {}", (converged) ? "TRUE" : "FALSE");
  long fcalls3 = matter_hop.getForceCalls();
  fcallsMove = fcalls2 - fcalls1;
  fcallsRelax = fcalls3 - fcalls2;
}
 
void GlobalOptimizationJob::randomMove(Matter &matter) {
  // create a random displacement
  AtomMatrix displacement;
  displacement.resize(matter.numberOfAtoms(), 3);
  displacement.setZero();
  int num = matter.numberOfAtoms();
 
  for (int i = 0; i < num; i++) {
    double disp = params.basin_hopping_options().displacement;
    if (!matter.getFixed(i)) {
      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");
          throw std::invalid_argument(std::format(
              "[Global Optimization] unknown displacement_distribution: {}",
              params.basin_hopping_options().displacement_distribution));
        }
      }
    }
  }
  matter.setPositions(matter.getPositions() + displacement);
}
 
void GlobalOptimizationJob::mdescape(Matter &matter) {
  int nmd;
  double ekinc, epot, etot, epot0, etot0;
  auto dyn = std::make_unique<Dynamics>(&matter, params);
  velopt(matter);
  epot = matter.getPotentialEnergy();
  ekinc = matter.getKineticEnergy();
  etot = ekinc + epot;
  epot0 = epot;
  etot0 = etot;
  size_t nummax = 0, nummin = 0;
  double enmin1 = 0.0, enmin2 = 0.0, en0000 = 0.0;
  double econs_max = -1.E100, econs_min = 1.E100, devcon;
  bool md_presumably_escaped = false;
  QUILL_LOG_DEBUG(log, "MD  {:5d}  {:20.10E}  {:15.5E}  {:15.5E}  ", 0,
                  epot - epot0, ekinc, etot - etot0);
  nmd = 1000;
  for (int imd = 1; imd <= nmd; imd++) {
    enmin2 = enmin1;
    enmin1 = en0000;
    dyn->velocityVerlet();
    epot = matter.getPotentialEnergy();
    ekinc = matter.getKineticEnergy();
    etot = ekinc + epot;
    en0000 = epot - epot0;
    if (enmin1 > enmin2 && enmin1 > en0000)
      nummax += 1;
    if (enmin1 < enmin2 && enmin1 < en0000)
      nummin += 1;
    QUILL_LOG_TRACE_L1(log,
                       "MD  {:5d}  {:15.5f}  {:15.5f}  {:12.2E}  {:4}  {:4}",
                       imd, epot - epot0, ekinc, etot - etot0, nummax, nummin);
    econs_max = std::max(econs_max, ekinc + epot);
    econs_min = std::min(econs_min, ekinc + epot);
    if (nummin >= static_cast<size_t>(mdmin)) {
      if (nummax != nummin)
        QUILL_LOG_WARNING(log, "WARNING: iproc,nummin,nummax {} {}", nummin,
                          nummax);
      md_presumably_escaped = true;
      break;
    }
  } // end of loop over imd
  devcon = econs_max - econs_min;
  if (md_presumably_escaped) {
    devcon = devcon / static_cast<double>(matter.numberOfFreeAtoms() * 3);
    if (devcon / ekin < 2.E-3) {
      ParametersLoadAccess::dynamics_options(params).time_step *= 1.1;
    } else {
      ParametersLoadAccess::dynamics_options(params).time_step /= 1.1;
    }
  } else {
    QUILL_LOG_DEBUG(log, "TOO MANY MD STEPS  ");
    ParametersLoadAccess::dynamics_options(params).time_step *= 2.0;
  }
}
 
void GlobalOptimizationJob::velopt(Matter &matter) {
  AtomMatrix vat(matter.numberOfAtoms(), 3);
  double tt1, tt2, tt3, vtot[3];
  int iat;
  vtot[0] = 0.0;
  vtot[1] = 0.0;
  vtot[2] = 0.0;
  for (iat = 0; iat < matter.numberOfAtoms(); iat++) {
    tt1 = eonc::rng::randomDouble();
    tt2 = eonc::rng::randomDouble();
    tt3 = eonc::rng::randomDouble();
    vat(iat, 0) = (tt1 - 0.5) * 2.0;
    vat(iat, 1) = (tt2 - 0.5) * 2.0;
    vat(iat, 2) = (tt3 - 0.5) * 2.0;
    vtot[0] += vat(iat, 0);
    vtot[1] += vat(iat, 1);
    vtot[2] += vat(iat, 2);
  }
  QUILL_LOG_DEBUG(log, "Linear momentum  {:15.5E}  {:15.5E}  {:15.5E}  ",
                  vtot[0], vtot[1], vtot[2]);
  vtot[0] /= matter.numberOfAtoms();
  vtot[1] /= matter.numberOfAtoms();
  vtot[2] /= matter.numberOfAtoms();
  for (iat = 0; iat < matter.numberOfAtoms(); iat++) {
    vat(iat, 0) -= vtot[0];
    vat(iat, 1) -= vtot[1];
    vat(iat, 2) -= vtot[2];
  }
  matter.setVelocities(vat);
  long nFreeCoords = matter.numberOfFreeAtoms() * 3;
  if (nFreeCoords <= 0) {
    throw std::invalid_argument("GlobalOptimizationJob::velopt: no free atoms");
  }
  double kinE = matter.getKineticEnergy();
  double kB = params.constants().kB;
  double kinT = (2.0 * kinE / nFreeCoords / kB);
  if (!(kinT > 0.0)) {
    throw std::runtime_error(
        "GlobalOptimizationJob::velopt: zero kinetic temperature");
  }
  double temperature = (2.0 * ekin / kB);
  matter.setVelocities(vat * std::sqrt(temperature / kinT));
}
 
void GlobalOptimizationJob::insert(Matter &matter) {
  double epot;
  size_t jlo, jlo_insert;
  epot = matter.getPotentialEnergy();
  jlo = hunt(epot);
  QUILL_LOG_DEBUG(log, "JLO= {}  {:10.5f}  ", jlo, std::abs(epot - earr[jlo]));
  if (!(std::abs(epot - earr[jlo]) <
        params.structure_comparison_options().energy_difference)) {
    jlo_insert = jlo;
    if (epot > earr[jlo])
      jlo_insert++;
    earr.insert(earr.begin() + jlo_insert, 1, epot);
  }
}
 
size_t GlobalOptimizationJob::hunt(double epot) {
  // epot is in interval [earr(jlo),earr(jlo+1)[ ; earr(0)=-Infinity ; earr(n+1)
  // = Infinity
  size_t jlo;
  double de;
  for (jlo = 0; jlo < earr.size(); jlo++)
    if (epot < earr[jlo])
      break;
  if (jlo == earr.size())
    jlo--;
  de = std::abs(epot - earr[jlo]);
  if (jlo > 0)
    if (std::abs(epot - earr[jlo - 1]) < de)
      jlo--;
  return jlo;
}
 
} // namespace eonc