ReplicaDynamicsJob.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/ReplicaDynamicsJob.h"
#include "eon/BaseStructures.h"
#include "eon/Dynamics.h"
#include "eon/ForceCallTimer.h"
#include "eon/HelperFunctions.h"
#include <stdexcept>
 
#include <format>
#include <fstream>
 
std::vector<std::string> ReplicaDynamicsJob::run() {
  auto seed = std::make_shared<Matter>(pot, params);
  std::string reactantFilename =
      eonc::helpers::getRelevantFile(params.main_options.conFilename);
  if (!eonc::io::io_ok(seed->con2matter(reactantFilename))) {
    QUILL_LOG_CRITICAL(log, "Failed to load {}", reactantFilename);
    throw std::runtime_error("failed to load " + reactantFilename);
  }
  (void)runFromMatter(std::move(seed));
  return returnFiles;
}
 
std::shared_ptr<Matter>
ReplicaDynamicsJob::runFromMatter(std::shared_ptr<Matter> initial) {
  if (!initial) {
    throw std::runtime_error("runFromMatter: initial Matter is null");
  }
  current = initial;
  current->setPotential(pot);
  reactant = std::make_shared<Matter>(pot, params);
  saddle = std::make_shared<Matter>(pot, params);
  product = std::make_shared<Matter>(pot, params);
  finalState = std::make_shared<Matter>(pot, params);
  finalStateTmp = std::make_shared<Matter>(pot, params);
 
  minimizeFCalls = mdFCalls = refineFCalls = dephaseFCalls = 0;
  time = 0.0;
 
  QUILL_LOG_DEBUG(log, "Minimizing initial reactant");
  {
    eonc::ForceCallTimer timer(minimizeFCalls);
    *reactant = *current;
    reactant->relax();
  }
 
  initExtra();
 
  int status = dynamics();
 
  saveData(status);
  reportResults();
 
  return current;
}
 
bool ReplicaDynamicsJob::checkState(Matter *curr, Matter *react) {
  Matter tmp(pot, params);
  tmp = *curr;
  tmp.relax(true);
  return !tmp.compare(*react);
}
 
long ReplicaDynamicsJob::refine(
    const std::vector<std::shared_ptr<Matter>> &buff, Matter *react) {
  QUILL_LOG_TRACE_L1(log, "Refining transition time.");
 
  long lo = 0;
  long hi = static_cast<long>(buff.size()) - 1;
 
  while ((hi - lo) > 1) {
    long mid = lo + (hi - lo) / 2;
    if (!checkState(buff[mid].get(), react)) {
      lo = mid;
    } else {
      hi = mid;
    }
  }
 
  return (lo + hi) / 2 + 1;
}
 
void ReplicaDynamicsJob::dephase() {
  long DephaseSteps =
      static_cast<long>(params.parallel_replica_options.dephase_time /
                        params.dynamics_options.time_step);
  Dynamics dephaseDynamics(current.get(), params);
  QUILL_LOG_DEBUG(log, "Dephasing for {:.2f} fs",
                  params.parallel_replica_options.dephase_time *
                      params.constants.timeUnit);
 
  long step = 0, loop = 0;
 
  while (step < DephaseSteps) {
    long dephaseBufferLength = DephaseSteps - step;
    loop++;
    std::vector<std::shared_ptr<Matter>> dephaseBuffer(dephaseBufferLength);
 
    for (long i = 0; i < dephaseBufferLength; i++) {
      dephaseBuffer[i] = std::make_shared<Matter>(pot, params);
      dephaseDynamics.oneStep();
      *dephaseBuffer[i] = *current;
    }
 
    bool transitionFlag = checkState(current.get(), reactant.get());
 
    if (transitionFlag) {
      long dephaseRefineStep = refine(dephaseBuffer, reactant.get());
      QUILL_LOG_DEBUG(log, "loop = {}; dephase refine step = {}", loop,
                      dephaseRefineStep);
      long ts = dephaseRefineStep - 1;
      ts = (ts > 0) ? ts : 0;
      QUILL_LOG_DEBUG(
          log,
          "Dephasing warning: in a new state, inverse the momentum and restart "
          "from step {}",
          step + ts);
      *current = *dephaseBuffer[ts];
      AtomMatrix velocity = current->getVelocities();
      velocity = velocity * (-1);
      current->setVelocities(velocity);
      step = step + ts;
    } else {
      step = step + dephaseBufferLength;
      QUILL_LOG_TRACE_L1(log, "Successful dephasing for {} steps", step);
    }
 
    if ((params.parallel_replica_options.dephase_loop_stop) &&
        (loop > params.parallel_replica_options.dephase_loop_max)) {
      QUILL_LOG_DEBUG(
          log,
          "Reach dephase loop maximum, stop dephasing! Dephased for {} steps",
          step);
      break;
    }
    QUILL_LOG_DEBUG(log, "Successfully Dephased for {:.2f} fs",
                    step * params.dynamics_options.time_step *
                        params.constants.timeUnit);
  }
}
 
void ReplicaDynamicsJob::saveData(int status) {
  std::string resultsFilename("results.dat");
  returnFiles.push_back(resultsFilename);
  size_t totalFCalls = minimizeFCalls + mdFCalls + dephaseFCalls + refineFCalls;
 
  {
    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("{} total_force_calls\n", totalFCalls);
      out << std::format("{} force_calls_dephase\n", dephaseFCalls);
      out << std::format("{} force_calls_dynamics\n", mdFCalls);
      out << std::format("{} force_calls_minimize\n", minimizeFCalls);
      out << std::format("{} force_calls_refine\n", refineFCalls);
      out << std::format("{} transition_found\n", (newStateFlag) ? 1 : 0);
 
      if (newStateFlag) {
        out << std::format("{:e} transition_time_s\n",
                           minCorrectedTime * 1.0e-15 *
                               params.constants.timeUnit);
        out << std::format("{:f} potential_energy_product\n",
                           product->getPotentialEnergy());
        out << std::format("{:f} moved_distance\n",
                           product->distanceTo(*reactant));
      }
 
      out << std::format("{:e} simulation_time_s\n",
                         time * 1.0e-15 * params.constants.timeUnit);
      out << std::format("{:f} speedup\n",
                         time / params.dynamics_options.steps /
                             params.dynamics_options.time_step);
    }
  }
 
  std::string reactantFilename("reactant.con");
  returnFiles.push_back(reactantFilename);
  if (!eonc::io::io_ok(reactant->matter2con(reactantFilename))) {
    QUILL_LOG_ERROR(log, "Failed to write {}", reactantFilename);
  }
 
  if (newStateFlag) {
    std::string productFilename("product.con");
    returnFiles.push_back(productFilename);
    if (!eonc::io::io_ok(product->matter2con(productFilename))) {
      QUILL_LOG_ERROR(log, "Failed to write {}", productFilename);
    }
 
    if (params.parallel_replica_options.refine_transition) {
      std::string saddleFilename("saddle.con");
      returnFiles.push_back(saddleFilename);
      if (!eonc::io::io_ok(saddle->matter2con(saddleFilename))) {
        QUILL_LOG_ERROR(log, "Failed to write {}", saddleFilename);
      }
    }
  }
}