eonc::GlobalOptimizationJob Class Reference

#include <GlobalOptimizationJob.h>

Inheritance diagram for eonc::GlobalOptimizationJob:

Public Member Functions
	GlobalOptimizationJob (std::unique_ptr< Parameters > parameters, Runtime &rt)
	~GlobalOptimizationJob (void)=default
void	hoppingStep (long, Matter &, Matter &)
void	decisionStep (Matter &, Matter &)
void	report (Matter &)
void	acceptRejectNPEW (Matter &, Matter &)
void	acceptRejectBoltzmann (Matter &, Matter &)
void	analyze (Matter &, Matter &)
void	examineEscape (Matter &, Matter &)
void	applyMoveFeedbackMD (void)
void	applyDecisionFeedback (void)
void	mdescape (Matter &)
void	randomMove (Matter &)
void	insert (Matter &)
size_t	hunt (double)
void	velopt (Matter &)
std::vector< std::string >	run (void)
	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.

Public Attributes
double	beta1
double	beta2
double	beta3
double	alpha1
double	alpha2
long	mdmin

Private Attributes
size_t	nlmin
long	fcallsMove
long	fcallsRelax
double	ediff
double	ekin
bool	firstStep
std::vector< double >	earr
std::string	escapeResult
std::string	decisionResult
std::string	hoppingResult
std::ofstream	monfile
std::ofstream	earrfile
eonc::log::Scoped	log

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 23 of file GlobalOptimizationJob.h.

Constructor & Destructor Documentation

GlobalOptimizationJob()

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

inline

Definition at line 25 of file GlobalOptimizationJob.h.

      : Job(std::move(parameters), rt),
        nlmin{0},
        ediff{1.E-1},
        ekin{5.E-2},
        beta1{params.global_optimization_options().beta},
        beta2{params.global_optimization_options().beta},
        beta3{1. / params.global_optimization_options().beta},
        alpha1{1. / params.global_optimization_options().alpha},
        alpha2{params.global_optimization_options().alpha},
        mdmin{params.global_optimization_options().mdmin},
        fcallsMove{0},
        firstStep{true},
        fcallsRelax{0},
        monfile{"monitoring.dat"},
        earrfile{"earr.dat"} {}

~GlobalOptimizationJob()

eonc::GlobalOptimizationJob::~GlobalOptimizationJob ( void )

default

Member Function Documentation

acceptRejectBoltzmann()

void eonc::GlobalOptimizationJob::acceptRejectBoltzmann	(	Matter &	matter_cur,
		Matter &	matter_hop )

Definition at line 219 of file GlobalOptimizationJob.cpp.

                                                                      {
  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";
  }
}

acceptRejectNPEW()

void eonc::GlobalOptimizationJob::acceptRejectNPEW	(	Matter &	matter_cur,
		Matter &	matter_hop )

Definition at line 209 of file GlobalOptimizationJob.cpp.

                                                                 {
  if (matter_hop.getPotentialEnergy() <
      matter_cur.getPotentialEnergy() + ediff) {
    decisionResult = "accepted";
  } else {
    decisionResult = "rejected";
  }
}

analyze()

void eonc::GlobalOptimizationJob::analyze	(	Matter &	matter_cur,
		Matter &	matter_hop )

Definition at line 76 of file GlobalOptimizationJob.cpp.

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

applyDecisionFeedback()

void eonc::GlobalOptimizationJob::applyDecisionFeedback

(

void

)

Definition at line 145 of file GlobalOptimizationJob.cpp.

                                                  {
  if (decisionResult == "accepted") {
    ediff *= alpha1;
  } else {
    ediff *= alpha2;
  }
}

applyMoveFeedbackMD()

void eonc::GlobalOptimizationJob::applyMoveFeedbackMD

(

void

)

Definition at line 124 of file GlobalOptimizationJob.cpp.

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

decisionStep()

void eonc::GlobalOptimizationJob::decisionStep	(	Matter &	matter_cur,
		Matter &	matter_hop )

Definition at line 185 of file GlobalOptimizationJob.cpp.

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

examineEscape()

void eonc::GlobalOptimizationJob::examineEscape	(	Matter &	matter_cur,
		Matter &	matter_hop )

Definition at line 111 of file GlobalOptimizationJob.cpp.

                                                              {
  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";
  }
}

hoppingStep()

void eonc::GlobalOptimizationJob::hoppingStep	(	long	istep,
		Matter &	matter_cur,
		Matter &	matter_hop )

Definition at line 244 of file GlobalOptimizationJob.cpp.

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

hunt()

size_t eonc::GlobalOptimizationJob::hunt

(

double

epot

)

Definition at line 411 of file GlobalOptimizationJob.cpp.

                                              {
  // 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;
}

insert()

void eonc::GlobalOptimizationJob::insert

(

Matter &

matter

)

Definition at line 396 of file GlobalOptimizationJob.cpp.

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

mdescape()

void eonc::GlobalOptimizationJob::mdescape

(

Matter &

matter

)

Definition at line 296 of file GlobalOptimizationJob.cpp.

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

randomMove()

void eonc::GlobalOptimizationJob::randomMove

(

Matter &

matter

)

Definition at line 266 of file GlobalOptimizationJob.cpp.

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

report()

void eonc::GlobalOptimizationJob::report

(

Matter &

matter_hop

)

Definition at line 153 of file GlobalOptimizationJob.cpp.

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

run()

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

virtual

Virtual run; used solely for dynamic dispatch.

Implements eonc::Job.

Definition at line 26 of file GlobalOptimizationJob.cpp.

                                                {
  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

velopt()

void eonc::GlobalOptimizationJob::velopt

(

Matter &

matter

)

Definition at line 352 of file GlobalOptimizationJob.cpp.

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

Member Data Documentation

alpha1

double eonc::GlobalOptimizationJob::alpha1

Definition at line 60 of file GlobalOptimizationJob.h.

alpha2

double eonc::GlobalOptimizationJob::alpha2

Definition at line 61 of file GlobalOptimizationJob.h.

beta1

double eonc::GlobalOptimizationJob::beta1

Definition at line 57 of file GlobalOptimizationJob.h.

beta2

double eonc::GlobalOptimizationJob::beta2

Definition at line 58 of file GlobalOptimizationJob.h.

beta3

double eonc::GlobalOptimizationJob::beta3

Definition at line 59 of file GlobalOptimizationJob.h.

decisionResult

std::string eonc::GlobalOptimizationJob::decisionResult

private

Definition at line 73 of file GlobalOptimizationJob.h.

earr

std::vector<double> eonc::GlobalOptimizationJob::earr

private

Definition at line 71 of file GlobalOptimizationJob.h.

earrfile

std::ofstream eonc::GlobalOptimizationJob::earrfile

private

Definition at line 76 of file GlobalOptimizationJob.h.

ediff

double eonc::GlobalOptimizationJob::ediff

private

Definition at line 68 of file GlobalOptimizationJob.h.

ekin

double eonc::GlobalOptimizationJob::ekin

private

Definition at line 69 of file GlobalOptimizationJob.h.

escapeResult

std::string eonc::GlobalOptimizationJob::escapeResult

private

Definition at line 72 of file GlobalOptimizationJob.h.

fcallsMove

long eonc::GlobalOptimizationJob::fcallsMove

private

Definition at line 66 of file GlobalOptimizationJob.h.

fcallsRelax

long eonc::GlobalOptimizationJob::fcallsRelax

private

Definition at line 67 of file GlobalOptimizationJob.h.

firstStep

bool eonc::GlobalOptimizationJob::firstStep

private

Definition at line 70 of file GlobalOptimizationJob.h.

hoppingResult

std::string eonc::GlobalOptimizationJob::hoppingResult

private

Definition at line 74 of file GlobalOptimizationJob.h.

log

eonc::log::Scoped eonc::GlobalOptimizationJob::log

private

Definition at line 77 of file GlobalOptimizationJob.h.

mdmin

long eonc::GlobalOptimizationJob::mdmin

Definition at line 62 of file GlobalOptimizationJob.h.

monfile

std::ofstream eonc::GlobalOptimizationJob::monfile

private

Definition at line 75 of file GlobalOptimizationJob.h.

nlmin

size_t eonc::GlobalOptimizationJob::nlmin

private

Definition at line 65 of file GlobalOptimizationJob.h.

---

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

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