eonc::ParallelReplicaJob Class Reference

#include <ParallelReplicaJob.h>

Inheritance diagram for eonc::ParallelReplicaJob:

Public Member Functions
	~ParallelReplicaJob ()=default
std::vector< std::string >	run () override
	Virtual run; used solely for dynamic dispatch.
std::shared_ptr< Matter >	runFromMatter (std::shared_ptr< Matter > initial)
	Matter-first entry; returns final trajectory Matter.
	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)
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.

Private Member Functions
void	dephase (Matter &trajectory, BondBoost *bias)
int	refineTransition (const std::vector< std::shared_ptr< Matter > > &snapshots, bool fake=false)

Private Attributes
std::vector< std::string >	returnFiles
std::shared_ptr< Matter >	reactant
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 ParallelReplicaJob.h.

Constructor & Destructor Documentation

~ParallelReplicaJob()

eonc::ParallelReplicaJob::~ParallelReplicaJob ( )

default

Member Function Documentation

dephase()

void eonc::ParallelReplicaJob::dephase ( Matter & trajectory, BondBoost * bias )

private

Definition at line 291 of file ParallelReplicaJob.cpp.

                                                                    {
  Dynamics dynamics(&trajectory, params);
  // The boost object names this Matter. Attach it here, and reset the
  // trajectory with assignKeepingBias so each dephase loop keeps it.
  trajectory.setBiasPotential(bias);
 
  const double dt = params.dynamics_options().time_step;
  if (!(dt > 0.0)) {
    throw std::invalid_argument(
        "ParallelReplicaJob::dephase: time_step must be positive");
  }
  int dephaseSteps = static_cast<int>(
      std::floor(params.parallel_replica_options().dephase_time / dt + 0.5));
  if (dephaseSteps < 1)
    dephaseSteps = 1;
  const long maxLoops =
      std::max(1L, params.parallel_replica_options().dephase_loop_max);
  QUILL_LOG_DEBUG(log, "[ParallelReplica] Dephasing: {} steps (max {} loops)",
                  dephaseSteps, maxLoops);
 
  Matter initial(pot, params);
  initial = trajectory;
 
  for (long loop = 0; loop < maxLoops; ++loop) {
    trajectory.assignKeepingBias(initial);
    dynamics.setThermalVelocity();
 
    for (int step = 1; step <= dephaseSteps; step++) {
      dynamics.oneStep(step);
    }
 
    Matter minimized(pot, params);
    minimized = trajectory;
    minimized.relax();
 
    if (minimized.compare(*reactant)) {
      QUILL_LOG_DEBUG(log, "[ParallelReplica] Dephasing successful");
      return;
    }
    QUILL_LOG_DEBUG(
        log,
        "[ParallelReplica] Transition occured during dephasing; Restarting");
  }
  // Exhausted retries: keep last dephased trajectory rather than hang.
  QUILL_LOG_DEBUG(log,
                  "[ParallelReplica] Dephase loop max reached; continuing");
}

Job() [1/4]

eonc::Job::Job ( std::shared_ptr< Potential > potPassed, const Parameters & parameters )

inline

Definition at line 91 of file Job.h.

      : jtype{parameters.main_options().job},
        params{parameters},
        owned_runtime_{std::make_unique<Runtime>()},
        runtime_{owned_runtime_.get()},
        pot{potPassed} {}

Job() [2/4]

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

inlineexplicit

Own a default-constructed Runtime.

Definition at line 89 of file Job.h.

      : Job(std::move(parameters), std::make_unique<Runtime>()) {}

Job() [3/4]

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

inline

Borrow: caller keeps Runtime alive (CLI stack / Python Session).

Definition at line 76 of file Job.h.

      : jtype{parameters->main_options().job},
        params{*std::move(parameters)},
        owned_runtime_{},
        runtime_{&rt},
        pot{helpers::makePotential(params.potential_options().potential, params,
                                   rt)} {}

Job() [4/4]

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

inline

Own a Runtime (one-shot makeJob / rvalue).

Definition at line 84 of file Job.h.

      : Job(std::move(parameters), *rt) {
    adoptRuntime(std::move(rt));
  }

refineTransition()

int eonc::ParallelReplicaJob::refineTransition ( const std::vector< std::shared_ptr< Matter > > & snapshots, bool fake = false )

private

Definition at line 339 of file ParallelReplicaJob.cpp.

                                                                  {
  int lo = 0;
  int hi = static_cast<int>(snapshots.size()) - 1;
 
  while ((hi - lo) > 1) {
    int mid = lo + (hi - lo) / 2;
    Matter snapshot(*snapshots[mid]);
    snapshot.relax(true);
 
    bool stillReactant;
    if (!fake) {
      stillReactant = snapshot.compare(*reactant);
    } else {
      stillReactant = static_cast<bool>(eonc::rng::randomInt(0, 1));
    }
 
    if (stillReactant) {
      lo = mid;
    } else {
      hi = mid;
    }
  }
 
  return (lo + hi) / 2 + 1;
}

run()

std::vector< std::string > eonc::ParallelReplicaJob::run ( )

overridevirtual

Virtual run; used solely for dynamic dispatch.

Implements eonc::Job.

Definition at line 28 of file ParallelReplicaJob.cpp.

                                             {
  reactant = std::make_shared<Matter>(pot, params);
  {
    const auto posIn =
        eonc::helpers::getRelevantFile(params.main_options().conFilename);
    if (!eonc::io::io_ok(reactant->con2matter(posIn))) {
      QUILL_LOG_CRITICAL(log, "Failed to load {}", posIn);
      throw std::runtime_error("failed to load " + posIn);
    }
  }
  (void)runFromMatter(reactant);
  return returnFiles;
}

runFromMatter()

std::shared_ptr< Matter > eonc::ParallelReplicaJob::runFromMatter

(

std::shared_ptr< Matter >

initial

)

Matter-first entry; returns final trajectory Matter.

Definition at line 43 of file ParallelReplicaJob.cpp.

                                                               {
  if (!initial) {
    throw std::runtime_error("ParallelReplicaJob::runFromMatter: null Matter");
  }
  reactant = initial;
  reactant->setPotential(pot);
 
  QUILL_LOG_DEBUG(log, "[ParallelReplica] Minimizing initial position");
  reactant->relax();
  if (!eonc::io::io_ok(reactant->matter2con("reactant.con"))) {
    QUILL_LOG_ERROR(log, "Failed to write reactant.con");
  }
 
  auto trajectory = std::make_shared<Matter>(pot, params);
  *trajectory = *reactant;
  Dynamics dynamics(trajectory.get(), params);
  BondBoost bondBoost(trajectory.get(), params);
 
  BondBoost *bias = nullptr;
  if (params.hyperdynamics_options().bias_potential ==
      Hyperdynamics::BOND_BOOST) {
    bondBoost.initialize();
    trajectory->setBiasPotential(&bondBoost);
    bias = &bondBoost;
  }
 
  dephase(*trajectory, bias);
 
  int stateCheckInterval = static_cast<int>(
      std::floor(params.parallel_replica_options().state_check_interval /
                     params.dynamics_options().time_step +
                 0.5));
  int recordInterval = static_cast<int>(
      std::floor(params.parallel_replica_options().record_interval /
                     params.dynamics_options().time_step +
                 0.5));
  if (stateCheckInterval < 1)
    stateCheckInterval = 1;
  if (recordInterval < 1)
    recordInterval = 1;
 
  std::vector<std::shared_ptr<Matter>> mdSnapshots;
  std::vector<double> mdTimes;
  double transitionTime = 0;
  Matter transitionStructure(pot, params);
  size_t refineForceCalls = 0;
 
  double simulationTime = 0.0;
  if (params.hyperdynamics_options().bias_potential == Hyperdynamics::NONE) {
    QUILL_LOG_DEBUG(
        log, "[ParallelReplica] {:>8} {:>12} {:>10} {:>12} {:>12} {:>10}",
        "Step", "Time (s)", "KE", "PE", "TE", "KinT");
  } else {
    QUILL_LOG_DEBUG(
        log,
        "[ParallelReplica] {:>8} {:>12} {:>10} {:>10} {:>12} {:>12} {:>10}",
        "Step", "Time (s)", "Boost", "KE", "PE", "TE", "KinT");
  }
 
  for (int step = 1; step <= params.dynamics_options().steps; step++) {
    if (params.hyperdynamics_options().bias_potential ==
        Hyperdynamics::BOND_BOOST) {
      // oneStep() evaluates accelerations (and therefore boost()) more than
      // once. Advance the equilibration counter here, once per MD step, so
      // ParallelReplica and SafeHyper share the same rmd_time schedule.
      bondBoost.advance();
    }
    dynamics.oneStep();
    double boost = 1.0;
    if (params.hyperdynamics_options().bias_potential ==
        Hyperdynamics::BOND_BOOST) {
      double boostPotential = bondBoost.boost();
      double kB = params.constants().kB;
      boost = std::exp(boostPotential / kB / params.main_options().temperature);
      simulationTime += params.dynamics_options().time_step * boost;
    } else {
      simulationTime += params.dynamics_options().time_step;
    }
 
    double kinE = trajectory->getKineticEnergy();
    double potE = trajectory->getPotentialEnergy();
    double kinT = (2.0 * kinE / (trajectory->numberOfFreeAtoms() * 3) /
                   params.constants().kB);
 
    if (step % params.debug_options().write_movies_interval == 0) {
      if (params.hyperdynamics_options().bias_potential ==
          Hyperdynamics::NONE) {
        QUILL_LOG_DEBUG(log,
                        "[ParallelReplica] {:>8} {:>12.4e} {:>10.4f} "
                        "{:>12.4f} {:>12.4f} {:>10.2f}",
                        step,
                        simulationTime * params.constants().timeUnit * 1e-15,
                        kinE, potE, kinE + potE, kinT);
      } else {
        double boostPotential = bondBoost.boost();
        QUILL_LOG_DEBUG(
            log,
            "[ParallelReplica] {:>8} {:>12.4e} {:>10.3e} "
            "{:>10.4f} {:>12.4f} {:>12.4f} {:>10.2f}",
            step, simulationTime * params.constants().timeUnit * 1e-15, boost,
            kinE, potE + boostPotential, kinE + potE + boostPotential, kinT);
      }
    }
 
    // Snapshots for refinement
    if (step % recordInterval == 0 &&
        params.parallel_replica_options().refine_transition) {
      auto snap = std::make_shared<Matter>(pot, params);
      *snap = *trajectory;
      mdSnapshots.push_back(std::move(snap));
      mdTimes.push_back(simulationTime);
    }
 
    // Check for transition
    if (step % stateCheckInterval == 0 ||
        step == params.dynamics_options().steps) {
      QUILL_LOG_DEBUG(log, "[ParallelReplica] Checking for transition");
 
      Matter minimized(pot, params);
      minimized = *trajectory;
      minimized.relax();
 
      if (!minimized.compare(*reactant) && transitionTime == 0) {
        QUILL_LOG_DEBUG(log, "[ParallelReplica] Transition occurred");
 
        if (params.parallel_replica_options().refine_transition &&
            !mdSnapshots.empty() && !mdTimes.empty()) {
          QUILL_LOG_DEBUG(log, "[ParallelReplica] Refining transition time");
          int snapshotIndex;
          {
            eonc::ForceCallTimer timer(refineForceCalls);
            snapshotIndex = refineTransition(mdSnapshots);
          }
          if (snapshotIndex < 0)
            snapshotIndex = 0;
          if (snapshotIndex >= static_cast<int>(mdSnapshots.size()))
            snapshotIndex = static_cast<int>(mdSnapshots.size()) - 1;
 
          transitionTime = mdTimes[static_cast<size_t>(snapshotIndex)];
          transitionStructure =
              *mdSnapshots[static_cast<size_t>(snapshotIndex)];
        } else {
          transitionStructure = *trajectory;
          transitionTime = simulationTime;
        }
        QUILL_LOG_DEBUG(log, "[ParallelReplica] Transition time: {:.3e} s",
                        transitionTime * params.constants().timeUnit * 1e-15);
        // A run that keeps going after the transition keeps the bond boost
        // it started with; plain copy assignment would clear it.
        trajectory->assignKeepingBias(transitionStructure);
        // A false stop_after_transition keeps the remaining dynamics steps.
        if (params.parallel_replica_options().auto_stop) {
          break;
        }
 
      } else if (step + 1 == params.dynamics_options().steps &&
                 transitionTime == 0) {
        // Fake refinement to prevent force-call bias (only if snapshots exist)
        if (params.parallel_replica_options().refine_transition &&
            !mdSnapshots.empty()) {
          QUILL_LOG_DEBUG(
              log,
              "[ParallelReplica] Simulation ended without seeing a transition");
          QUILL_LOG_DEBUG(
              log, "[ParallelReplica] Refining anyways to prevent bias...");
          {
            eonc::ForceCallTimer timer(refineForceCalls);
            refineTransition(mdSnapshots, true);
          }
        }
        transitionStructure = *trajectory;
      }
 
      mdSnapshots.clear();
      mdTimes.clear();
    }
  }
 
  std::unique_ptr<Matter> product;
  if (transitionTime != 0) {
    int decorrelationSteps = static_cast<int>(
        std::floor(params.parallel_replica_options().corr_time /
                       params.dynamics_options().time_step +
                   0.5));
    if (decorrelationSteps < 0) {
      decorrelationSteps = 0;
    }
    QUILL_LOG_DEBUG(log, "[ParallelReplica] Decorrelating: {} steps",
                    decorrelationSteps);
    for (int dstep = 1; dstep <= decorrelationSteps; ++dstep) {
      dynamics.oneStep(dstep);
    }
    product = std::make_unique<Matter>(pot, params);
    *product = *trajectory;
    product->relax();
    if (!eonc::io::io_ok(product->matter2con("product.con"))) {
      QUILL_LOG_ERROR(log, "Failed to write product.con");
    }
  }
 
  // Write results
  std::string resultsFilename("results.dat");
  returnFiles.push_back(resultsFilename);
  {
    std::ofstream out(resultsFilename, std::ios::binary);
    if (out) {
      out << std::format(
          "{} potential_type\n",
          magic_enum::enum_name<PotType>(params.potential_options().potential));
      out << std::format("{} random_seed\n", params.main_options().randomSeed);
      out << std::format("{:f} potential_energy_reactant\n",
                         reactant->getPotentialEnergy());
      out << std::format("{} force_calls_refine\n", refineForceCalls);
      out << std::format("{} total_force_calls\n",
                         PotRegistry::get().total_force_calls());
 
      if (transitionTime == 0) {
        out << "0 transition_found\n";
        out << std::format("{:e} simulation_time_s\n",
                           simulationTime * params.constants().timeUnit *
                               1.0e-15);
      } else {
        out << "1 transition_found\n";
        out << std::format("{:e} transition_time_s\n",
                           transitionTime * params.constants().timeUnit *
                               1.0e-15);
        out << std::format("{:e} correlation_time_s\n",
                           params.parallel_replica_options().corr_time *
                               params.constants().timeUnit * 1.0e-15);
        out << std::format("{:f} potential_energy_product\n",
                           product->getPotentialEnergy());
      }
      out << std::format("{:f} speedup\n",
                         simulationTime /
                             (params.dynamics_options().steps *
                              params.dynamics_options().time_step));
    }
  }
 
  // bondBoost is a stack local of this function and setBiasPotential does not
  // own it, so the trajectory must not carry the pointer past the return.
  // getBiasForces reads a null bias potential as zero bias, which is what a
  // Matter that is no longer being boosted means.
  trajectory->setBiasPotential(nullptr);
 
  return trajectory;
}

Member Data Documentation

log

eonc::log::Scoped eonc::ParallelReplicaJob::log

private

Definition at line 38 of file ParallelReplicaJob.h.

reactant

std::shared_ptr<Matter> eonc::ParallelReplicaJob::reactant

private

Definition at line 33 of file ParallelReplicaJob.h.

returnFiles

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

private

Definition at line 32 of file ParallelReplicaJob.h.

---

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

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